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ABSTRACT 

We study the relation between the density distribution of tracers for large-scale struc- 
ture and the underlying matter distribution - commonly termed bias - in the ACDM 
framework. In particular, we examine the validity of the local model of biasing at 
quadratic order in the matter density. This model is characterized by parameters b\ 
and 62. Using an ensemble of ./V-body simulations, we apply several statistical meth- 
ods to estimate the parameters. We measure halo and matter fluctuations smoothed 
on various scales. We find that, whilst the fits are reasonably good, the parameters 
vary with smoothing scale. We argue that, for real-space measurements, owing to the 
mixing of wavemodes, no smoothing scale can be found for which the parameters are 
independent of smoothing. However, this is not the case in Fourier space. We measure 
halo and halo-mass power spectra and from these construct estimates of the effective 
large-scale bias as a guide for b\ . We measure the configuration dependence of the halo 
bispectra Bhhh and reduced bispectra Qhhh for very large-scale fc-space triangles. From 
this data we constrain b\ and &2, taking into account the full bispectrum covariance 
matrix. Using the lowest-order perturbation theory, we find that for i?hhh the best-fit 
parameters are in reasonable agreement with one another as the triangle scale is var- 
ied; although, the fits become poor as smaller scales are included. The same is true for 
Qhhh- The best-fit values were found to depend on the discreteness correction. This led 
us to consider halo-mass cross-bispectra. The results from these statistics supported 
our earlier findings. We then developed a test to explore whether the inconsistency in 
the recovered bias parameters could be attributed to missing higher-order corrections 
in the models. We prove that low-order expansions are not sufficiently accurate to 
model the data, even on scales k\ ~ 0.04 /iMpc -1 . If robust inferences concerning 
bias are to be drawn from future galaxy surveys, then accurate models for the full 
nonlinear bispectrum and trispectrum will be essential. 
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1 INTRODUCTION 

The accurate estimation and modelling of higher-order 
clustering statistics in current and future galaxy redshift 
surveys has the potential to act as a powerful probe 
for cosmological physics. The higher-order connected mo- 
ments, beginning at lowest order with the three-point 
correlation function and its Fourier analogue, the bis- 
pectrum, when interpreted within the gravitational insta- 
bility paradigm, encode important information regarding 
the growth of large-sca le structure (|Matarrese et al.lll997l ; 
IScoccimarro et all 1993). Their measurements also provide 
insight into the statistical nature of the primordial fluc- 
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tuations jFrv fc Scherrerlll994l; ISefusatti fc Komatsul [20071 : 



iNishimichi et al.ll2010l: iBaldauf et alj|201lf ) and the cosmo- 
logical parameters |Sefusatti et al.ll2006l h~Another attribute 
of three-point statistics, and the focus of our study, is their 
capability to probe the manner in which an observable 
tracer population of objects, such as galaxies, is related to 
the unobservable matter distribu t ion - termed the 'bias ' 
(iKaiserl Il984l: iDekel fc Reesl Il987l; |Frv fc Gaztanagal Il993l : 



iDekel fc Lahavll 19991 : ICatelan et alJl200Ch 



If the primo rdial fluctuations were Gaussian, as appears 
to be the case jKomatsu et alj|2010t l. then the statistical 
properties of the initial fields are fully characterized by the 
power spectrum, with all higher-order connected correlators 
vanishing. However, gravitational instability leads to the 
coupling of Fourier modes and this generates a hierarchy 
of non-vanishing connected correlators, each of which has 
a precise characteristic mathematical structure. The matter 
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bispectrum is thus an inherently nonlinear quantity, whose 
signal depends on closed triangles in Fourier space. In the- 
ory, this should vanish at early times and on scales large 
enough where linear theory is valid. If galaxy bias is local 
and linear, then the bispectrum of the observable tracers 
is proportional to the matter bispectrum. If on the other 
hand, bias is local and nonlinear then the triangle configu- 
ration dependence of the signal is modified, and this happens 
in a very precise and calculable w ay. Thus the bispectrum 
can be used to constrain the bias ((Fry fe Gaztanagalll993l ; 
iMatarrese et al.lll997l ; IScoccimarro et al.lll999l '). 

There is a long and rich history of measurements of 
three-point s tatistics from ga l axy su rveys, going all the 
way back to iPeebles fc Grothl (1 19751 ). However, attempts 
to constrain the nonlinearity of galaxy bias from galaxy 
redshift surveys have only been performed over the last 
decade. iFeldman et all (|200lh and lScoccimarro et all (|200ll ) 
both analyzed the IRAS survey using the bispectrum, and 
found a negative quadratic bias; although, d ue to small 
sampl e size the constraints were rather weak. IVerde et alj 
(2002) analyzed the 2dFGRS survey, also using the bispec- 
trum approach, and claimed that the flux-limited sample 
was an unbiased tracer of the dark matt er. A subsequent 
analys is of the final 2dFGRS data set by iGaztanaga et all 
|2005h . using the 3-point correlation function, contradicted 
this: using information from weakly non-linear scales (R ~ 
6 — 27 h~ 1 Mpc) the unbiased case (6i = 1 and b 2 = 0) 
was excluded at the order of 9a. More recently a num- 
ber of authors have analyzed various data releases of the 
SDSS (iNishimichi et al.l|2007l - DR3. iMcBride et all [ioill 
- DR6 and iMarinl |201o[ - DR7). These all claim a non- 



zero quadratic bias term for most of the samples within the 
dataset. Obviously, these variations of the results with sur- 
vey and statistical method require an explanation. 

Whilst the local bias model can be used to test whether 
the bias is linear or nonlinear, a significant detection of non- 
zero nonlinear bias does not imply that we have understood 
the bias. In order to believe that these measurements are 
meaningful, we need to be sure that the local model is indeed 
the correct model for interpreting data. This is currently an 
open question. Attempting to shed light on this subject is 
one of the aims of this paper. Over the past few years, the 



of authors (Heavens et al. 


19981; Gaztanaga & Scoccimarro 


2005 


; ISmith et al.1 2007; 


Guo & Jing 20091; IManera et al. 


201C 


; iRoth & Porcianil l201ll; iManera & Gaztanagal l201ll). 



However, no firm conclusions have yet been reached. 

In this paper, we use a large ensemble of 40 mid- 
resolution, large volume, pure dark matter iV-body simu- 
lations, to test the validity of the local bias model. In this 
study we compare a selection of different methods for deter- 
mining the bias. We first present a point-wise comparison 
of the halo and matter density fields smoothed on certain 
scales. We also utilize the power spectrum to estimate a 
large-scale effective bias parameter. Then we expend most 
of our efforts on using the bispectrum and reduced bispec- 
trum approach for constraining the bias. Besides the auto- 
bispectra, we also present, for the first time, measurements 
of the halo-matter cross-bispectra: Bhhm and Bt mm , and 
Qhhm and Qhmm- The value of these new statistics becomes 
apparent when correcting for shot-noise effects. Finally, we 
perform a numerical test that allows us to sharply illumi- 



nate the importance of terms in the theory that are beyond 
the tree-level expansions typically used. 

The paper is divided up as follows: in 32]we present the 
theory for the matter bispectrum and the local bias model. 
In 33 we provide details of the numerical simulations used 
in this work. In 3H we present the results for the bias pa- 
rameters from the various commonly used simple estimators. 
In 32] we present the estimation of the bias from the bispec- 
trum. Then, in 3H]we present measurements of bias from the 
cross-bispectra. In 33 we present a test of the importance of 
terms in the theory that are beyond tree level. Finally in 3H 
and 3H] we discuss our findings and present our conclusions. 



2 THEORETICAL OVERVIEW 

2.1 Standard perturbation theory dynamics 

In the fluid approximation, the gravitational collapse of col- 
lisionless cold matter structures in the expanding Universe, 
can be fully characterized by specifying the evolution of the 
density 5p(x) and the p eculiar velocity <5v(x) perturbations 
(|Bernardeau et ai1l2002T ). Focusing primarily on the density 
field, we work with models of the matter density contrast: 

p(x,t) -p(t) 



P(t) 



(1) 



where p(t) is the mean matter density of the Universe. In 
Fourier space, we define its corresponding Fourier represen- 
tation, <5(k), accordingly, as 



5(x 



k r <5(k)e-*' x 5(k) 



d 3 xc5(x)e ikx .(2) 



(2tt) 3 

It can be shown that the nonlinear equations of motion 
for 8(k) can be solved exactly by perturbative expansions 
of th e type l|juszkiewiczlll98 li : IVishniad[l983l : iGoroff et ail 
I1986T ): 



5(k) = (t)5„(k) , 

n = l 

and, where <5 n (k) is given by, 



*n(k) = 



d 3 qi 

(2tt) 3 
xF n (qi, 



d 3 q„ 
(2tt) 3 
,q»)<5i(qi 



(3) 

(27r) 3 ,S D (k-q 1 q n ) 

*i(q«) • (4) 

The density kernel F n is the dimensionless, homoge- 
neous, mode coupling function that couples together the 
amplitudes and phases of n initial Fourier wavemodes 
{f)(qi), . . . , <5(a n )}. As was shown by dGoroff et al.l 1 19861 : 
iMakino et ail 1 19921 ; I Jain fc Bertschingerj ll994T l. the nth 
order kernel may be constructed recursively from the 
lower-order solutions. Linear theory is thus represented by 
•Fi(qi) = 1, and the first nonlinear correction by F%, where 



v n 1 -\ 5 ki ■ k 2 
F 2 (ki,k 2 ) = - + 

7 K1K2 



fci fe 
fe fe 



2 (ki-k 2 )' 

7 /c-j^ k*2 



(5) 



The above approach defines the standard perturbation the- 
ory (hereafter SPT). Before moving on, we note that the 
above statements are only exactly true for the Einstein- 
de Sitter model. However, it has been shown that the 
F2 kernel is almost independent of c osmology (jFrvl 1 19941 ; 
iBouchet et~al]|l995l ; rHivon et al.lll995l ). We therefore adopt 
Eq. ((5} when dealing with the density at second order. 
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2.2 From dynamics to statistics 

Owing to the stochastic nature of the density field, we are 
not interested in reproducing a specific density field per se, 
but instead in characterizing its statistical properties. In this 
work we focus on 2- and 3-point correlation functions in 
Fourier space. These we may write as: 

(tf(ki)J(ka)) = (2^) 3 5 D (k 12 )P mm (k 1 ) ; (6) 
(<5(ki)5(k 2 )<5(k 3 )) = (27r) 3 5 D (k 123 )B mmm (k 1 ,k 2 ,k 3 ), (7) 

where P mm (k) and B mmm (ki, k 2 , k 3 ) constitute definitions 
of the power and bispectrum. For the Dirac delta functions 
we used the short-hand notation <5 D (ki.., n ) = 5 I3 (ki + . . . + 
k n ) and these guarantee that P and B are translationally in- 
variant. This is an important property for estimation, since 
it means that we should consider only closed pairs and tri- 
angles in Fourier space: ~}2 ki = 0. 

The perturbative expansion of the density field de- 
scribed in the previous section implies that P mm and B mmm 
may also be described in a perturbative fashion. Hence, 

(*(ki)*(ka)) = ([<5i(ki) + < 5 2 (k 1 ) + ...] 

x [<5i(k 2 ) + 5 2 (k 2 ) + ...]) ; (8) 
(<5(ki)...<5(k 3 )) = ([«5i (ki) + J 2 (ki) + ...]... 

x [5i(k 3 )+<5 2 (k 3 ) + ...]) . (9) 

Since we are assuming that the initial Fourier modes are 
Gaussianly distributed, i.e. the phase of each initial mode 
is uniformly random <j> £ [0, 2ir], modes must cancel in 
pairs. Hence, Wick's theorem applies, and so odd prod- 
ucts of initial Fourier modes must vanish: (5i(ki)iS 2 (k 2 )) = 
(<5i(ki)<5i(k 2 )<5i(k 3 )) = . This leads us to write the per- 
turbative expansions for P mm and -Bmmm as: 



Pmm 

B 

m mm 

(ki,k 2 ) 



— Pnm(fc) + Pmm(^) + • ■ • 

= B^ m (ki,k 2 ) + BL 1 im(ki,k 2 ) 



(10) 



(11) 



We shall refer to the lowest order terms in the expansions as 
'tree-level' terms. For P, P' ' is simply the linear spectrum, 
while for B the tree-level term can be written: 

S m l m (ki,k 2 ) = 2Pmi(A ; i)P m i(fe)P2(k 1 ,k 2 ) + 2cyc .(12) 

In this work we shall mainly be dealing with tree-level quan- 



and B 



tities; we now set: Pmm = Px 
otherwise stated. 

Another statistical quantity commonly used to 
explore galaxy clustering is the reduced bispectrum 
|Peebles fe Groth|[l975l : IScoccimarro et al.ll 19981 ). which can 
be defined: 



Qmmm(kl, k 2 , k 3 



Pmmm(kl, k 2 , k 3 ) 
ni(fa)Pmm(fe) + 2 CyC 



(13) 



As will be made clear below, the importance of this statistic 
becomes apparent when one considers non-Gaussian terms 
that are generated by simple quadratic products of Gaussian 
fields. In this case Q mmm simply scales as a constant. 



only formed in dark matter haloes 
tion for all models of galaxy formation (|White 



as is the us ual assumr . 

ReeJl97a 



then understanding the clustering of haloes is an essential 
comp onent of any theory of galaxy biasing (|Smith et al.l 
120071 ). In the local model of halo biasing, the number den- 
sity of dark matter haloes of mass scale M, smoothed over 
a scale R, can be expressed as a function of the local matter 
density, smoothed on the sc ale R. This function may then 
be Taylor expanded to give dFrv fc Gaztanagal 1 19931 : 1 Coles! 
Il993l : IMP et al.|[l997l : ISmith et al.ll2007l ): 



<5 h (x|M, R) = Y, 
3=0 



[*(x|ft)] 



(14) 



where we defined the smoothed halo over-density to be, 
6h(x\M,R) = [n h (x|M, R) -n h (M)]/n h (M). Owing to the 
fact that (5h) = 0, the constant coefficien t bo(M) = 
- Y,%2 b i( M ) ( &j ) /i ! (|Frv fc Gaztanagal Il993l ). Note that 
on Fourier transforming <5h(x|M, R) the constant 60 only 
contributes to the unmeasurable k — mode. The terms 
61 (M) and 6 2 (M) represent the linear and first nonlinear 
bias parameters, respectively. 

In Fourier space, Eq. (|14p can be written as: 



S h (k\M,R) = bi(M)5(k\R) 



+ 



b 2 (M) f d 3 qi 
(2*) 



r <5(qi|P)<5(k-qi|P) + 



,(15) 



where Si(qj\R) = W(\qj\R)5i(qj). If one inserts the SPT 
expansions for the density into the local model, then, up to 
second order in the density and bias, one finds: 

5 h (k\M,R) = fci(M) [<5i(k|P) + <5 2 (k|P)] 

+ ^/05 1 (q 1 |^ 1 (k-q 1 ]P).(16) 

Using this approach one may then find a perturbative ex- 
pansion for the halo power and bispectra: 



P hh (M) = P^(M) + P^(M) + 



Bhhh(M) = p<° h ' h (M) + B,^ h (M) + 



(i) 



(17) 
(18) 



Again, we refer to the lowest order terms in these expansions 
as tree-level terms, and for these we have: 



r^(k\M) 



b\{M)Vmm{k) 



(19) 



B 



(0) 



hh (ki,k 2 |M) = 6?(M)B mmm (ki,k2) + 6i(M)6 2 (M) 

X [p mm (fcl)Pmm(fc 2 ) + 2 CyC ] , (20) 



where in the above expressions we have derived the spec- 
tra of the smoothed fields: V = W 2 (kR)P(k), and B = 
W (kiR)W (k2R)W (ksR) B . However, when we estimate the 
bispectrum from data we do not smooth the fields apart 
from the CIC assignment scheme used to obtain the density 
cont rast field . As p ointed out by ISmith et all (|2007i . 12008 1 
and ISefusattl |2009l ), one way to overcome this is to adopt 
the ansatz: 



2.3 Halo Bias: Local Form 

In this study we investigate the relation between the clus- 
tering of dark matter haloes and total matter. If galaxies are 



B™ h (ki,k 2 ,k 3 |M) = 



W 2 (kR) ' 

(ki,k 2 ,k 3 |M,i?) 
W(k 1 R)W(k 2 R)W(k 3 R) 



(21) 
(22) 



© 0000 RAS, MNRAS 000, 000-000 



4 J. E. Pollack, R. E. Smith & C. Porciani 



On applying this 'de-smoothing' operation to Eqs (|19|l 
and (|20|l . one finds: 



P^(k\M) = fei(M)P mm (fc) ; 



B^ki.kalM) 



(23) 

6?(M)B mmm (ki,k 2 ) + 6 1 (M) 6 2 (Af) 

,k2 ^ram (&1 )-Pmm(^2 ) 4" 2 CyC . 

(24) 



where we have defined the function (jSefusa tti 2009): 
W(\\n\R)W(\k 2 \R) 



Hi 



ki,k 2 



W(|kl + k 2 |i?) 



(25) 



Note that in the limit of very large scales or arbitrarily small 
smoothing scales, kiR — > for i € {1, 2, 3}, Eq. Q24p approx- 
imates to: 



£C h (ki,k 2 |M) 



bl (M)B mmm (ki , ka) + 6? (Af ) &a (AO 
x [P mm (fci)P mm (fc 2 ) + 2 eye ] . (26) 



Again, since in this paper we are only considering tree- level 
expressions we shall take p£?- — > Phh and B^L — > -Bhhh- 

Considering now the reduced halo bispectrum, it may 
be defined in a similar fashion to Eq. (|13p : 



Qhhh(ki,k 2 ,k 3 |Af) 



B h hh(ki,k 2 |M) 



Phh(fci|M)P hh (fc 2 |M)+2cyc 



(27) 



On inserting our tree-level expressions from Eqs ()23|l 
and (l24l , we find that 



Qhhh(M) = 



Q 



&i(M) b\{M) 
where we have defined, 



. 62 (Af) . , . 
+ , 771 fN a(ki, k 2 ,k 3 ) , 



a(ki, k 2 , k 3 



Wk 1 ,k 2 P(fci)P(fe) + 2cyc 



(28) 



(29) 



P(fei)P(fc 2 ) + 2cyc 

Again, in the limit of very large scales or arbitrarily small 
smoothing scales and a(ki,k 2 ,k3) — > 1, the above expres- 
sion approximates to: 

mmm(kl, k 2 , k 3 ) 6 2 (Af ) 



Q hhh (ki,k 2 ,k 3 |Af) 



bi(Af) 



6?(Af) 



(30) 



We now see the utility of the reduced bispectrum: if one con- 
structs halo/galaxy density fields from a local transforma- 
tions of the matter density, then the lowest order nonlinear 
corrections will lead to a function that is a scaled version of 
the matter Qmmm, plus a constant offset. Moreover, if the 
density field were simply Gaussian, then estimates of Qhhh 
on large scales would directly measure fe 2 /fei. 



3 iV-BODY SIMULATIONS 

For our investigations of the bias, we use an ensemble of 
40 large iV-body simulations, executed on the zBOX-2 and 
zBOX-3 supercomputers at the University of Zurich. We use 
only the z — outputs from the simulations. Each simula- 
tion was performed using the publicly available Gadget-2 
code (|Springelll20"05l ), and followed the nonlinear evolution 
under gravity of N = 750 3 equal-mass particles in a comov- 
ing cube of length L s im = 1500/i _1 Mpc. 

The cosmological model that we simulated was analo- 
gous to the basic va nilla ACDM model de termined by the 
WMAP experiment (|Komatsu et al.ll2009l ): matter density 



fi m = 0.25, vacuum density Qa = 0.75, power spectrum nor- 
malization ag = 0.8, power spectral index n = 1, and dimen- 
sionless Hubble parameter h — 0.7. The transfer function 
for the simulations was generated using the publicly avail- 
able cmbf ast code jSeliak fc Zaldarriaga|[l996l : ISeliak et al.l 
I2003T I. with high sampling of the spatial frequencies on large 
scales. Initial conditions were set at redshift z — 49 us- 
i ng the serial vers i on of the publicly available 2LPT code 
(|Scoccimarrdll998l : ICrocce et al.ll2006l 1. 

Dark matter halo catalogues were generated for each 
simulation using the Friends-of-Friends (FoF) algorithm 
(|Davis et al.lll985h . with the linking- length parameter b — 
0.2, where b is the fraction of the inter-particle spacing. For 
this we employed the fast parallel B-FoF code, provided to 
us by V. Springel. The minimum number of particles for 
which an object is considered to be a bound halo was set 
at 20 particles. This gave a minimum host halo mass of 
Afmin = 1.11 x 10 13 /i -1 Me. For our analysis of the bias, 
we use the full sample of haloes and this corresponded to 
roughly N h ~ 1.26 x 10 6 haloes per sim ulation. Furth er de- 
tails of the simulations may be found in ISmithl |2009h . 



4 SIMPLE ESTIMATES OF BIAS 

Before we examine halo bias in the context of the bispec- 
trum, we explore two alternative methods for studying the 
bias. We first evaluate the second-order local biasing model 
directly, by comparing in a point-wise fashion the halo and 
matter density fields, smoothed over a range of scales. Then, 
we use power spectra to determine an effective large-scale 
bias. 



4.1 Analyzing Density Fields 

One obvious way to examine the local model of biasing is 
to simply construct a scatter plot of the local density of 
dark matter haloes against the local density of dark matter 
in the simulations (s ee for example ISheth fc Lemsonl Il999l : 
iDekel fc Lahavlll999l ). As was discussed in i]2.3l this model 
only makes sense in the context of smooth fields. We shall 
therefore also inspect how the model parameters depend on 
the adopted smoothing scale R. 

We generate the smoothed density fields as follows: we 
assign particles/haloes to a Fourier grid using the CIC algo- 
rithm (c.f. 35TT}; then we Fourier transform the grid using 
the FFT algorithm; each Fourier mode is then multiplied by 
a Gaussian filter of the form: 



W{kR) = exp \-(kR) 2 /2] 



(31) 



Finally, on taking the inverse Fourier transform, we obtain 
the smoothed S(x\R) and S h (x\R). We perform the above 
procedure for 28 of the ensemble of simulations and consider 
the filter scales: R = {50 , 20 , 10} /i _1 Mpc. 

In Figure [1] we present the bin averaged scatter plots 
of <5 h (x|f?) vs. <5(x|f?), averaged over the realizations. The 
colour contours are shaded by the normalized population 
density of that pixel, e.g. the central red region indicates that 
most of the points in the simulation are regions of density 
close to average. We also see that as the smoothing scale is 
decreased (panels going from left-to-right), that the scatter 
increases and that there are more points that have higher 
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Figure 1. Scatter plots of <5h(x) versus <5(x) smoothed with a Gaussian filter of various scales averaged over the realizations. From left 
to right, the panels correspond to the smoothing scales R = {50, 20, 10} h~ 1 Mpc. The color coding denotes the log of the population 
density, i.e. the red region corresponds to the largest concentration of points and the white background to null values. The dot-dashed 
line in each panel denotes the local halo bias model up to second-order with the best-fitting bias parameters averaged over 28 realizations. 



R 


bo ± -"b 


fel ± cr bl 


-2 ± o"b 2 


[h^Mpc] 


xl0~ 3 






50 


1.3 ± 0.1 


1.497 ± 0.002 


-0.577 ± 0.031 


20 


12.0 ± 0.1 


1.542 ± 0.006 


-0.635 ± 0.004 


10 


37.2 ± 0.1 


1.644 ± 0.005 


-0.512 ± 0.001 



Table 1. Average of the mean bias parameters and the root-mean 
square errors for the local halo bias model up to second-order 
averaged over 28 realizations determined from fitting the scatter 
plots of the halo and matter density fields smoothed on scales 
k s = {0.02, 0.05, 0.1} /iMpc" 1 . 



and lower density. Conversely, as the filter scale is increased 
the relation becomes tighter and more linear. One obvious 
conclusion that may be drawn from this behaviour is that 
the bias relation is certainly not deterministic. 

In order to obtain a more quantitative understanding, 
we next consider fitting for the parameters of the local bias 
model at second-order. From Eq. (|14p we have: 

.h(x|M, R) = b (M) + b 1 (M) 6{x\R) + ^p- {5{x\R)} 2 .(32) 

We perform a least-squares analysis on each realization, and 
then average over the resulting set of bias parameters to ob- 
tain the mean parameters: bo(M), fei(M), and &2(M). The 
lcr-errors are then estimated in the usual way, as quadratic 
deviations from the sample mean. In Fig. [T] we plot the re- 
sultant best-fit local model as the dot-dashed line in each of 
the three panels. 

The information on the parameters is summarized in 
Table [T] as a function of the filter scale R. This clearly 
shows that the estimates of b\ increase as the smooth- 
ing scale is decreased, whereas those for 62 appear to be 
parabolical. Naively, one might expect that the nonlinear 
bias terms should approach zero as the amount of smooth- 
ing is increased and nonlinearities are washed out, however, 
at R = 50/i _1 Mpc even with a(x\R) < 1, the fluctuations 



are still significant enough to yield a non-zero -2- Note also 
that in all cases 60 7^ 0. 

The local model, as written in Eq. (|14[) . asserts that the 
parameters bi are independent of the smoothing scale R, and 
we, therefore, consider the implications as follows. Suppose 
that nonlinear bias is exactly as described by Eq. (|32[) . but 
that the coefficients are not independent of the smoothing 
scale. Let us now consider the results that would be obtained 
from measurements for two smoothing scales R a and R b . 
From Eq. (|32[) we would have: 

5 h (x|i? ) = &S + b?5(x|i4) + y [-(x|i? Q )] 2 ; (33) 

5 h (x\R b ) = 6f, + & 1 5(x|i? i) ) + ^[<5(x|%)] 2 . (34) 

Supposing now that we desmoothed each of the fields, by 
Fourier transforming and dividing out the appropriate win- 
dow function. We would then have: 

5 h (k) = 6^(k) + ||-^5(q)5(|k-q|)W q , k - q (i?a) 

5 h (k) = 6^(k) + ||-^5(q)5(|k-q|)lV q , k - q (J? 6 ) . 

In order for the above equations to be equivalent, then we 
must have 

K = b\ (35) 

h a h b \ W(qRb) W{\k-Cj\R a )W{kR a ) } 

2 2 [W(qRa) W(\k-q\R b ) W(kR b )\ ' ( ' 

The last of the two equations may only be satisfied if and 
only iiR a = R b or {kR, qR, \k-q\R} < 1. Since the <5 h (x|7?) 
vs. (5(x|i?) method is inherently a real space measure it in- 
volves contributions from all Fourier modes. It is therefore 
difficult to ensure that b% = b\. 

We conclude that the above method will not be a safe 
way to recover bias parameters independent of the smooth- 
ing scale. We now turn to Fourier space methods. 
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Figure 2. Scale dependence of effective bias parameters b^h' 6Sr and fe^ m (c.f. Eqs (138 I t and ]39ll ). estimated from the auto- and 

cross-power spectrum as a function of wavemode. For the left and central panels: solid blue and open red symbols denote the bias when 
Phh is not and is shot noise corrected, respectively. The first panel shows b^h when the nonlinear matter power spectrum is used; the 
second panel shows the same but when the linear matter power spectrum is used; the third panel shows bhmi where the red stars and 
blue points denote the case where the nonlinear and linear matter power spectra are used, respectively. 



fcJ/iMpr 1 ] 




,NL,SC 
°hh 






0.03-0.09 


1.589 ± 0.002 


1.493 ± 0.002 


1.589 ± 0.004 


1.487 ± 0.004 


0.04-0.12 


1.624 ± 0.002 


1.486 ± 0.002 


1.638 ± 0.003 


1.489 ± 0.003 


0.05-0.15 


1.663 ± 0.002 


1.474 ± 0.002 


1.709 ± 0.003 


1.503 ± 0.003 


0.06-0.18 


1.695 ± 0.001 


1.460 ± 0.002 


1.775 ± 0.002 


1.511 ± 0.003 



Table 2. Weighted average estimates of the effective bias, ft^h (see Eq. (138 I t for a defintions). 



4.2 Effective large-scale bias from power spectra 

We now use various halo power spectra to derive estimates 
for an effective large-scale halo bias. 

In order to do this, we first measure the Fourier trans- 
form of the matter and halo density fields as described in 
AppendixfA] The halo-halo, halo- mass and mass-mass power 
spectra, {Phh, -Phm, fmm}, are then estimated from the data 
by performing the following sums: 

v w(*0 

h»{h) = j^t £ SMStQst) , (37) 

where {fj,, v\ £ {m, h}, is the sample volume (which in 
our case is the simulation volume), and where N(k) are the 
number of Fo urier modes i n a sh ell of thickness Afc. 

Following [Smith et all (|2007m . we next construct the es- 
timators: 




where N s is the number of simulations and is the linear 
matter power spectrum. Note that in the case of bhh we also 
consider shot-noise corrected versions of these two estima- 
tors, i.e. we correct Phh using Eq. (|Allfl . We denote these 

NL SC L SC 

bias estimates by b hh ' and b hh , respectively. Finally, we 
determine the la errors by evaluating the variance of each 
realization against the mean. 



The first panel of Fig. [2] shows b^t (solid blue data 
points) and bj^' 30 (open red points). The bias for the shot- 
noise corrected terms remains roughly constant at ~ 1.49 
down to scales k ~ 0.08 h Mpc -1 and with very small er- 
rors, indicating that the result is highly constrained by the 
data. On scales smaller than this the bias is a decreasing 
function of k. Without shot-noise correction, we find that 
the bias is strongly scale dependent, and the bias rapidly 
increases with increasing k. 

The second panel of Fig.[2] shows the results obtained 
from using (solid blue points) and fe|j h SC (open hexago- 
nal points). The results are similar to those for b^h"' 80 , but 
with increased cosmic variance on large scales. An oscillation 
st ructure is also prese nt, this can be understood as explained 
in iGuzik et al.l(|2007l l. Nevertheless, comparing the two pro- 
vides a clear indication of the validity of the tree-level power 
spectrum up to k = 0.08 ftMpc -1 . 

The third panel of Fig. [2] shows the bias results fr^m 
(solid red points) and fo^ (solid blue points). The value 
of b^m stays roughly constant for the whole scale range 
considered in the estimate, while b^ m is not smooth and 
clearly shows the imprint of the oscillation structure. Nev- 
ertheless, we find bh m ~ 1.48, to within the errors for 
k < 0.08 /iMpc -1 . On comparing the results for b[^' sc 
and bj^, we see that for scales k < 0.08 feMpc -1 , these esti- 
mates are compatible and that the effective large-scale bias 
is roughly b ~ 1.49. Interestingly, these findings are consis- 
tent with real space measures of the effectiv e large scale bias 
from cell variances (jSmith fc Marianll201lh . 

In Table [2] we report the weighted average and cor- 
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responding la error on the effective bias, &hh, computed 
over the same fc-modes corresponding to the magnitude of 
the third wavevector kz for each triangle configuration. The 
tabulated results for the analysis of the uncorrected data 
confirms the results shown in Figure [2] that bias is indeed 
scale-dependent. Applying the shot-noise correction yields a 
more constant effective bias, even for the range of fc-modes 
entering into our bispectrum estimation. Interestingly, the 
value 61 = 1.49 ± 0.002 found for k G [0.03, 0.09], is in good 
agreement with the result for 61 obtained from fitting the 
density fields smoothed on scales R = 50/i _1 Mpc. There- 
fore, if we opt to infer that the effective bias is equivalent 
to 61 over these scales, then the bispectrum (reduced bis- 
pectrum) should also yield this value for 61 when fitted over 
the same scale ranges (c.f. Tabled}. That is, if the local bias 
model is correct and the tree-level bispectrum (reduced bis- 
pectrum) is a sufficient description of the nonlinearities on 
these scales. 

Before moving on, we point out that one can also use 
the h alo power spectra to define an effective 62 (|Smith et alj 
l2009h . however we shall not explore this here. 



5 HALO BIAS FROM BISPECTRA 

In this section, we present our main results from the analysis 
of the halo bispectra. 



5.1 Bispectrum Estimation 

The computational code used to estimate the matter and 
halo power a nd bispectra i s a m odified version of the code 
developed in ISmith et al.l (120081). which itself is based on 
the algorithm of Scoccimarro et all l|l998l ). The major mod- 
ification to that code, which we have implemented, is that 
no random subsampling of the Fourier modes is performed 
to estimate the bispectrum. Instead, all modes that con- 
tribute to a particular triangle configuration are used. In 
this work we use a FFT grid of size N e = 512 3 to esti- 
mate the power and bispectra. We only evaluate triangles 
that have = 2fci, but consider the variation of B with 
the angular separation of the two vectors. The largest scale 
at which we estimate the bispectrum is fci = 0.03 /iMpc -1 , 
and this is » 7.5kt, where fef = 2n/L « 0.004 h Mpc -1 . Fur- 
ther details of the bispectrum estimation procedure may be 
found in Appendix 1X1 

Figure [3] shows the ensemble-averaged shot-noise cor- 
rected results for the halo bispectra £?hhh (open red squares) 
and matter bispectra B mmm (solid blue diamonds), mea- 
sured from the ensemble of iV-body simulations. The four 
panels show the results obtained for the scales: k\ = 
{0.03,0.04,0.05, 0.06} /iMpc -1 . The error bars are the 1-a 
errors on the mean, derived from the ensemble to ensemble 
variation. The solid red line represents the tree-level pre- 
diction for Bmmm as given by Eg. (fl2)l . We see that this 
appears to be a good description of the -Bmmm estimates for 
the scales that we have considered. We notice that, for the 
case fci = 0.06 /iMpc -1 , the theory systematically under- 
predicts the measurements for 9 /it ~ 0.5 (but see ^5.41 for a 
more quantitative discussion of the goodness of fit). 

Figure [4] shows the same as in Fig.[3] however this time 
for Qhhh and Q 

mmm • 



5.2 Bias estimation 

Our method for estimating the bias p arameters foll o ws an 
app roach similar to that present ed by IScoccimarr 3 (120001 ) 
and H orciam fc Giavaliscol (|2002i l. To start, we take a \ 
function that is a quadratic form of the type: 

Ng Ng 

X 2 (6i,6 2 ) =J2J2 Ai ( bl ^y^ A ^ bl ^ > ( 4 °) 

i=l 3=1 

where Ne is the number of angular bins considered and 



A s = 



-Bhhh(fci, k2, 0i) — gghh (fejj k2,0i\bi, 62) 



(41) 



Chhh(fel, k2,0i) 

Note that in dividing the difference between the estimate of 
the ensemble average bispectrum (B) and the model predic- 
tion (-B mod ) by the standard deviation (o"hhh), » — 1 ij is in fact 
the inverse correlation matrix. Recall that the correlation 
and covariance matrices are related by: ry = Cy / \JCaCjj . 

In order to minimize this x 2 function and so recover the 
best-fit bias parameters, we need an estimate of 1 — 1 y, and 
we do this using the standard unbiased estimator: 



1 ^Bto 

lm l > fc =i 



13 ~ (N si 

where TVsim is the number of simulations and where 
^ (k) l (*i,kA)-Eu*(*i>*a»0i) 



A _ ^hhhV 



<Thhh(fei, fea, 9i 



(42) 



(43) 



where in the above is the k tn estimate of the bispec- 

trum and where Bhhh is the mean of the ensemble. 

Next, we use singular-value-decomposition (SVD) to in- 
vert the rij. For the estimate of the inverse correlation ma- 
trix, we utilize principal component analysis (PCA) to re- 
move some of the noisy eigenvectors. We select the fraction 
of principal components that account for 95% of the vari- 
ance. According to this selection criteria we typically retain 
15 out of 20 of the most 'dominant' eigenmodes. Note that 
as pointed out in lHartlap et al.1 (|2007f ). C" 1 / How- 
ever, since we are using PCA, this should be a subdominant 
correction. Thus, we may approximate Eq. (|40[) as: 

N g Ng 

X 2 = ^^A,(6 1 ,6 2 )[ii T A J R]r. 1 A ;J (6 1 ,62) , 

i=i 3=1 

Ng Ng Ng 

= ^^A i (6i,6 2 )^^A^ 1 iZ y A,(6i,6 2 ) , 

i=l 3=1 1 = 1 

Ng 

« ^Ki'YfOu, (44) 
;=i 

where the correlation matrix r was diagonalized by rota- 
tion into its eigenbasis, i.e. r = R T AR, with A repre- 
senting a diagonal matrix of eigenvalues. We also defined 
Yi = X^j=i RuAi(bi, 62)- Note that in the final approximate 
expression we include a matrix On, this a diagonal matrix 
with entries either 1 or 0, depending on whether the eigen- 
vector is to be retained or cut from the PCA reconstruction. 

Finally, the X 2 {bi, fe) function was minimized using the 
Levenberg-Marquardt routine for non-linear least squares 
fitting. 



th 
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Figure 3. Ensemble-averaged matter and halo bispectrum measurements for 40 LCDM TV-body simulations in real-space in comparison 
with the PT models at tree-level. Each panel shows the shot-noise corrected bispectrum measurements as a function of angle for a variety 
of triangle configurations at different scale ranges set by fci = {0.03 0.04 0.05 0.06} ftMpc -1 and k 2 = 2fci. The blue solid symbols 
represent the matter bispectrum. whereas the open squares denote the halo bispectrum. The tree-level bispectrum is represented by the 
solid orange line and the local halo bias model with the best- fitting parameters listed in Table HI is the dashed violet line. 



5.3 Errors in parameter estimates 

To the best-fit parameters (6i, 62), we assign both systematic 
and statistical errors. 

In our context, the systematic errors correspond to the 
errors induced in the best-fit parameters from fitting the 
data with a noisy inverse covariance matrix (or correlation 
matrix). Owing to the relatively low number of simulations 
(A^im = 40), we expect that Eq. (|42|) provides a noisy es- 
timate of rZ ■ In order to estimate the errors this has on 
the best-fit parameters we employ the jackkni fe subsampling 
method (see for example lNorberg et al.|[2009l ). This involves 
slicing the total data set into iV su b subsamples. Then a re- 
sampling of the data is obtained by excluding one of the sub- 
samples from the set. From this resampling we then estimate 
the mean statistic of interest and the inverse correlation ma- 
trix as described in the previous section. The resampled data 
set is then used to determine a new estimate of the best-fit 
bias parameters. This procedure is then repeated for all of 



the possible 7V su b resamplings of the data. In our particu- 
lar case we treat the measurements from each simulation as 
the regions to be included or excluded, and this gives us 40 
jackknife estimates of the bias parameters (61,62). The pa- 
rameter covar iance matrix for the systematic errors can be 
computed as jNorberg et alj|200g| ): 



N sut 



N sub 



bj)(b. 



■bi) , 



(45) 



where 6;^ is the estimate of bi from the fc th resampling of 
the data and bi is the estimate of the mean bi obtained from 
all of the resamplings. 

The statistical error is obtained directly from the non- 
linear least-squares analysis. The routine mrqmin provides 
an approximation to the errors on the parameters that cor- 
responds to a Ax 2 ~ 1 for a one-parameter model. However, 
the confidence regions we present in the forthcoming plots 
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Figure 4. Ensemble-averaged matter and halo reduced bispectrum measurements of the 40 LCDM ./V-body simulations in real-space in 
comparison with the PT models at tree-level. Point and line styles as in Figure [3] 



correspond to either Av 2 = (2.30,6.17), which roughly de- 
note the (~ lcr, ~ 2a) errors for a two-parameter model. 

Given that we consider the two forms of error: system- 
atic and statistical), and that one is never consistently larger 
than the other, in all forthcoming tables, we choose to re- 
port only the total error. This is obtained simply from the 
two errors added in quadrature. 



5.4 Testing the validity of the tree-level matter B 

Before we report the estimates of the halo bias parameters, 
we first present a test of the validity of the tree-level model 
for the matter bispectrum. We do this by applying the y 2 
test described above, to the B mmm and Q mmm data, and 
so fit for 6i and &2- Note that since the total number of 
principal components retained equals 15, then for a two- 
parameter model the number of degrees-of-freedom equals 
13. If the tree-level expressions in the large-scale limit as 
given by Eqs (|26p and (|30p are correct, then we should expect 
to find b\ — 1 and c*2 = 0. 

Table [3] presents the best-fit nonlinear bias parameters 



for the four different bispectrum scale ranges discussed ear- 
lier. In the analysis we fit the shot-noise corrected bispectra. 
The y 2 values (end column of the table), confirms that the 
tree- level expressions Bmmm an( i Qmmm! provide good fits for 
the triangle configurations with ki = {0.03, 0.04} h Mpc -1 . 
However, for ki = {0.05, 0.06} h Mpc" 1 the fits are poor 
given the \ 2 estimates, and we see that, for both B and Q, 
they yield non-zero values for &2 at lcr. The results also im- 
ply that the failure of the tree-level model on these scales 
is more severe for Q than for B. This can be understood 
by noting that b\ from Q mmm shows a prominent departure 
from unity, whereas Bmmm does not (although the deviation 
still exceeds 2a). We thus conclude that it is likely that the 
tree-level expressions for the halo bispectra will only be valid 
for fci sC 0.04/i Mpc" 1 , for our chosen bispectra configura- 
tions. 



5.5 Constraints on b\ and 62 from halo bispectra 

Table [4] presents the best-fit nonlinear bias parameters and 
their respective la errors in quadrature, obtained from the 
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ki [hMpc- 1 ] 




bi ± o- i)1 


b 2 ± f(, 2 


x 2 


0.03 


Ommm 
Qnunm 


1.01 ± 0.07 
0.93 ± 0.19 


-0.04 ± 0.25 
-0.05 ± 0.30 


19.08 
19.03 


0.04 


Bmnmi 
Qmmm 


0.98 ± 0.03 
1.05 ± 0.09 


0.04 ± 0.10 
0.07 ± 0.16 


14.31 
14.95 


0.05 


Qmmm 


0.97 ± 0.02 
1.14 ± 0.07 


0.13 ± 0.08 
0.19 ± 0.13 


38.83 
19.97 


0.06 


Qmmm 


0.98 ± 0.02 
1.15 ± 0.04 


0.10 ± 0.10 
0.22 ± 0.08 


34.09 
29.58 



Table 3. Assessment of the validity of the tree-level modelling 
by fitting the matter bispectra and reduced bispectra. Column 1: 
bispectrum triangle scale; column 2: statistic, where -B mmm and 
Qmmm are shot-noise corrected; columns 3 and 4: best-fit b\ and 
b 2 along with lrx errors; column 5: \ 2 - 



X 2 analysis of £>hhh and Qhhh- Note that we present the 
results for both the uncorrected and shot-noise corrected 
measurements, indicated in the table by superscript 'SC. 
Table [4] also shows the \ 2 value of these best-fit parameters 
as an indication of the goodness-of-fit. 

In Figures [3] and [4] we also show the tree-level theo- 
retical models for Bhhh and Qhhh (dashed lines), where the 
best-fit bias parameters from Table[4]have been used. These 
figures demonstrate that, at least by-eye, the tree-level mod- 
els provide a reasonable description of the data. However a 
more detailed inspection of Table [4] reveals some important 
discrepancies. 

For the case of fitting B, the shot-noise correction is less 
important, as we see that the estimates of b\ for all bispectra 
configurations with and without shot-noise corrections are 
consistent to within the errors, and have b\ ~ 1.4. However, 
b 2 shows systematic differences, being more negative if the 
correction is made, and for this we find that b 2 ~ —0.25. On 
the other hand, for the case of Q, the results clearly show 
that the shot-noise subtraction has an important effect on 
the recovered values for the bias parameters. If the shot- 
noise is not corrected, then we see that the estimates for 61 
increase systematically as we go from triangle configurations 
with ki = 0.03 ftMpc" 1 to fei = 0.06ftMpc -1 . Whereas if 
it is corrected, then we find 61 ~ 1.8 and 62 ~ —0.3 to 
within the errors. On comparing the results from B and Q 
we see that, whilst the values for bi disagree significantly, 
surprisingly those for 62 remain consistent at the la level. 

The \ 2 function of Eq. (|41[) may be interpreted as 
a Gaussian likelihood if we make the transformation, 
£({.Bhhh}|&i,&2) oc exp[— x 2 /2]. Once suitably normalized 
and on assuming a set of prior probabilities we may then 
explore the shape of the confidence regions in the posterior 
probability p(bi , b 2 \ {B hhh }) . 

Figure [5] shows the la likelihood confidence contours 
in the posterior probability for the nonlinear bias parame- 
ters for the four scales considered according to our method 
of analysis described above. The solid lines denote the size 
of the confidence regions at the 68% and 95% level (i.e. 
Ax 2 ~ 2.3,6.17) when we construct a correlation matrix 
from the 40 realizations without regard for the systematic 
uncertainty. The dashed-lines demonstrate the magnitude 



fci [riMpc" 1 ] 




61 ± cr bl 


b 2 ± cr b2 


x 2 


0.03 


-Bhhh 
rSC 

Qhhh 
^ hhh 


1.43 ± 0.11 
1 A o 1 nil 

2.09 ± 0.55 
1.75 ± 0.47 


-0.18 ± 0.40 
-U.oO ± U.00 
-0.12 ± 0.76 
-0.39 ± 0.56 


17.20 
1 / .Uo 
16.08 
16.66 


0.04 


-Bhhh 
RSC 
-Bhhh 

Qhhh 
Vhhh 


1.41 ± 0.08 
l.oo ± U.Uo 
2.32 ± 0.39 
1.80 ± 0.29 


-0.05 ± 0.26 
-U.z ( ± U.zO 
0.14 ± 0.58 
-0.34 ± 0.36 


26.92 

Of? 1 P. 

ZO. 10 
26.40 
26.96 


0.05 


-Bhhh 
R SC 
-"hhh 

Qhhh 
hhh 


1.40 ± 0.06 
1.00 m u.uo 
2.66 ± 0.26 
1.90 ± 0.19 


0.15 ± 0.21 
-u.zo nz u.±o 
0.57 ± 0.42 
-0.30 ± 0.22 


31.92 

IZ.Do 

11.53 

11.60 


0.06 


-Bhhh 

B sc 

-°hhh 
Qhhh 
^hhh 


1.41 ± 0.05 
1.37 ± 0.03 
2.84 ± 0.20 
1.87 ± 0.14 


0.19 ± 0.24 
-0.23 ± 0.13 
0.88 ± 0.40 
-0.30 ± 0.19 


63.64 
19.30 
20.70 
19.47 



Table 4. Best-fit bias parameters from fitting the halo-halo-halo 
bispectra and reduced bispectra. Column 1: bispectrum triangle 
scale; column 2: statistic, where -Bhhh an d Qhhh are raw, and 

columns 3 and 4: 

„2 



where B^f h and 



^hhh are shot-noise corrected; 
best-fit bi and b 2 along with la errors; column 5: 



A " 



at which the 68% and 95% confidence regions expand fol- 
lowing our generation of a set of jackknife subsamples to 
monitor the effect due to the implicit error associated with 
the estimated correlation matrix. Hence, we clearly see the 
relevance of accounting for the uncertainty of the correlation 
matrix when obtaining the bias parameter constraints. The 
discrepancy between the resulting jackknife error ellipses for 
B and Q is less severe than the likelihood contours obtained 
from the complete sample where the level of agreement im- 
proves progressing to large scales, yet this might be due to 
the fact that the statistical error is more prominent at larger 
scales. Interestingly, the overlap of the two likelihood regions 
at 2a for fci=0.03, 0.04 and 0.05 /iMpc" 1 occurs with the 
rectangular region or strip denoting the effective bias mea- 
sure, &£J ' sc , at la. These set of panels in Figure [5] convey 
pictorially the information obtained from the results in Ta- 
ble [4] that the likelihood contours from analysis of Q show 
an evolution with decreasing scale toward larger and larger 
61, whereas the constraints on b 2 remain consistent. Lastly, 
the constraints obtained from analyzing the Q-amplitudes 
are much weaker than those coming from the bispectrum. 



6 HALO BIAS FROM CROSS-BISPECTRA 



In £|5 . 51 we saw that shot-noise corrections influenced the 
recovery of the bias parameters particularly for Q. In this 
section, we attempt to develop the use of cross-bispectra, as 
measures of the bias that are less susceptible to discreteness 
effects. The use of cross-correlations in large-scale structure 
work has lo ng been known as a way of reducing discreteness 
corrections ((Peebles 1980). However, it is only relatively re- 
cent t hat it has been applied to study bias rtSmith et al 
20071 ; iDalal et all 120081 ; ISmithl [20091; IPadmanabhan et al 



200g| ; lDesiacaues et al.ll2009l ; IPillepich et al.ll2ofoT ) 
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Figure 5. Evolution of the likelihood contours for the bias parameters bi and b 2 , estimated from -Bhhh and Qhhhi with scale. The 
solid lines denote the 68% and 95% confidence intervals, obtained from a full exploration of the likelihood surface around the best 
fit values; the dashed lines denote the same, but where the jackknife parameter covariance matrix from Eq. 1145 \ has been used to 
determine the error contours. The top left, top right, bottom left and bottom right panels show the results for triangle configurations 
with fci = {0.03, 0.04, 0.05, 0.06} /iMpc -1 , respectively. The vertical black lines denote the effective bias parameter &j^' SC , using the 
same wavemodes that enter into the bispectrum estimates. 



6.1 Definitions and theory 

We may define the halo cross-bispectra as follows: 

<<5 h (ki)«5 h (k 2 )5(k 3 ))= (27r) 3 <5 D (k 12 3)B hhm (ki,k 2 ,k 3 ) ;(46) 
(*h(ki)J(k 2 )«5(k 3 ))= (2tt) 3 <5 d (ki 23 )B h mm(ki,k2,k3) (47) 

We then symmeterize these quantities by the operations: 

^hhm^ = [-Bhhm + -Bhmh + B m hh] /3 ; (48) 
= [Bhmm + Bmhm + fimmh] /3 . (49) 

For ease of notation we shall now simply take B^™' = Bhhm 

and Bfo^m = Bhmm, unless otherwise indicated. We may 
now also define the cross-reduced bispectra as: 

Qhhm = Bhhm/PPhhm ', (50) 
Qhmm — B\i mm / P P\ lmm , (51) 



where we have for the denominators (again symmetrized): 
PPbhm = | [P hh (fci)P hm (fc 2 ) + 2 eye] 

+ i [ftm(fe)flm(fe) + 2 eye] ; (52) 

PPhmm = | [Phm (kl ) Pmm (fc 2 ) + 2 CyC ] 

+ i[P hm (fci)P hm (A: 2 ) + 2cyc] . (53) 

The relations for PPhhm and PPhmm can easily be con- 
structed using a graphical approach. Let us consider three 
nodes two of which are the same and the third is different 
(we shall think of the nodes as the density fields). Label these 
nodes 1, 2, and 3. Now consider all possible ways to connect 
the three nodes together by two edges. When two nodes, 
which are the same, connect together this gives us an auto- 
power spectrum with a delta function, and when two nodes 
that are different connect together this gives us a cross- 
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power spectrum and delta function. One may then sym- 
metrize the results by considering all possible relabellings 
of the nodes and dividing by three. 

In Appendix [B] we calculate the tree-level cross- 
bispectra, Bhmm and Bhhm, in the local model of halo bi- 
asing. The main results are: 

B h ?L(ki,k 2 ,k 3 ) « &i(M)BWU(ki,k 2 ,k 3 ) 
62 (M) 



[w kl 



k 2 J 



>(o) 



(kx)P^l(k 2 ) + 2 eye 



(54) 



B^^.ka.ka) 



6?(M)B m 'mm(k 1 ,k 2 ,k 3 ) + ~&i(M) 



xb 2 (M) [w kl , k2 Pi? m (fa)^i° m (fe) + 2cyc] . (55) 

In the limit of large scales, and or small smoothing scales, 
the filter functions M / k 1 ,k 2 ~~ ► 1 an d we have 

s£L « kB£L + - [P m ] m (fc 1 )P m , m (fc 2 ) + 2cyc] ; (56) 

^hhm « «ML» + -&1&2 [p4 t m(fcl)P4°m(fc2) + 2 CyC ] . 

(57) 

At second order in the nonlinear bias and PT, the cross- 
reduced bispectra are: 

mmm 6 2 (M) 



)(0) 
- hmm 



Q (0) 



2 + 61 (M) + 26i (M) + &?(M) 



a(ki,k 2 ,k 3 ) ;(58) 



+ 



2& 2 (M) 



26i (M) + 1 2&f (M) + 61 (M) 



a(ki,k 2 ,k 3 ) 



(59) 



In the large-scale limit, a — I 1, these expressions become 

o (0) _ 3Qmmm , b 2 (M) 



fhmm ~ 2 + fci(M) + 26i(M) + 6?(M) ' 

•-i(O) ^ 3 Qmmm 

26 2 (M) 

Ahm ~ 2&i(M) + 1 26? (M) + 6i(M) 



(60) 
(61) 



6.2 Estimation of the cross-bispectra 

The cross-bispectra Bhhm and Bh mm can be estimated fol- 
lowing the algorithm described in i]5.1l with some small mod- 
ifications. Firstly, the estimates must be symmetrized, and 
for the discrete form of Bhhm we have, 

1 V 2 Nt " 

B^ m (ki, fa, 0i2 ) = frf- Y] 

3 iVtri f— ' 

x I fte[5h(k ni )5h(k„ 2 )5m(k n3 )] + 2cyc| ,(62) 

and with a similar relation for Bh mm . The reduced bispectra 
are estimated by dividing the above bispectrum estimates 
by estimates for PPhhm an d PPhmm from Eqs Q52p and (1531) , 
respectively. 

One further complication is constructing the corrections 
for shot-noise. This may be p er forme d following the counts- 
in-ce lls approach of I Peebles! (|l980t ) (and see also ISmithl 
120091 ). We find that the symmetrized corrections for Bhhm 
and Bhmm can be written as: 



By 



3^ [ P ^ (k 



+ 2 eye 



Bhmm.shot = [Phm(fa) +2cyc| , 

in m L 



(63) 
(64) 



fci [hMpc- 1 ] 






&2 i °"b2 


x 2 




^hmra 


1.37 ± 1.18 


-0.36 ± 0.32 


16.86 


0.03 


dSC 
hmm 


Lot ± l.lo 


-U.oD ± U.oZ 


1 P. Q(Z 
10. OO 


li m 111 


1.98 ± 0.75 


-0.75 ± 0.99 


19.09 




Qf c 

" hmm 


1.98 ± 0.75 


-0.75 ± 0.99 


19.09 




-^hmm 


1.33 ± 0.16 


-0.33 ± 0.51 


16.83 


0.04 


dSC 
hm m 


l.oo ± U.10 


-U.00 ± U.ol 


ID. OO 


Qlinim 


2.50 ± 0.56 


-0.31 ± 0.89 


17.96 




Qf c 

c li m 111 


2.51 ± 0.56 


-0.32 ± 0.89 


17.97 




D 

^hmm 


1.30 ± 0.07 


-0.01 ± 0.23 


26.98 


0.05 


dSC 
hmm 


i.ou m u.u* 


n 09 -i- o 91 

-U.UZ Hi U.ZO 




Qhmm 


2.94 ± 0.25 


0.24 ± 0.57 


16.90 




Qf c 

" hmm 


2.94 ± 0.25 


0.23 ± 0.57 


16.95 




D 

-'-'hmm 


1.29 ± 0.07 


-0.001 ± 0.28 


13.99 


0.06 


-B sc 

hmm 


1.29 ± 0.07 


-0.004 ± 0.28 


13.93 


Qhram 


3.45 ± 0.30 


1.43 ± 0.94 


46.60 




Qf c 

^ hmm 


3.45 ± 0.30 


1.42 ± 0.94 


46.87 



Table 5. Best-fit bias parameters from halo-mass-mass bispec- 
tra and reduced bispectra. Column 1: bispectrum triangle scale; 
column 2: statistic, where -Bhmm and Qhmm are raw and where 
and Q^mm are shot-noise corrected; columns 3 and 4: best- 



BSC 

hmm 

fit 61 and 62 along with l<x errors; column 5: \ 2 



ki [hMpc- 1 ] 




61 ± a bl 


f> 2 ± <Jb 2 


x 2 


0.03 


-Bhhm 

nnm 
Qhhxa 

o sc 


1.42 ± 0.25 
1.42 ± 0.24 
2.08 ± 0.40 
2.11 ± 0.40 


-0.32 ± 0.60 
-0.45 ± 0.47 
-0.49 ± 0.46 
-0.67 ± 0.42 


17.13 
17.49 
17.45 
17.46 


0.04 


-Bhhm 
nnm 
Qhhxa 

o sc 


1.39 ± 0.14 
1.37 ± 0.15 
2.49 ± 0.46 
2.56 ± 0.49 


-0.24 ± 0.25 
-0.37 ± 0.24 
-0.17 ± 0.56 
-0.47 ± 0.51 


21.29 
21.42 
22.92 
24.60 


0.05 


-Bhhm 

rsc 

-°hhm 
Qhhm 

osc 


1.38 ± 0.08 
1.36 ± 0.07 
2.87 ± 0.18 
2.92 ± 0.19 


-0.05 ± 0.18 
-0.29 ± 0.14 
0.22 ± 0.29 
-0.37 ± 0.25 


22.46 
13.16 
12.55 
16.77 


0.06 


-Bhhm 

rsc 

-"hhm 
Qhhxa 

o sc 

^lihm 


1.39 ± 0.07 
1.36 ± 0.07 
3.12 ± 0.27 
3.25 ± 0.31 


-0.07 ± 0.25 
-0.31 ± 0.22 
0.54 ± 0.52 
-0.14 ± 0.47 


20.56 
15.32 
37.81 
60.54 



Table 6. Best-fit bias parameters from halo-halo-mass bispectra 
and reduced bispectra. Column 1: bispectrum triangle scale; col- 
umn 2: statistic, where -Bhhm and Qhhm are raw and where BSC m 

columns 3 and 4: best-fit b\ 



and Qj^hm are s hot-noise corrected; 
and fe 2 along with la errors; column 5: x 2 



where n m = N/V^ and nh = Nh/Vp, are the number den- 
sity of matter particles and haloes, respectively. For the re- 
duced bispectra we must correct the estimates of PPhhm and 
PPhmm, which are written in the following form: 

2 
3n h 
2 

3n n 

As it is the case that n m 2> rTh, we expect that the shot- 
noise corrections to Bhmm will be significantly smaller than 



Qdcnom 
hhm.shot 



-jdenom 
/hmm. shot 



3=^ [Phm(fci) + 2cyc 
[Phm(fci) + 2cyc] 



(65) 
(66) 
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Figure 6. Evolution of the 95% likelihood contours for b\ and hi obtained from the halo auto- and cross-bispectra and reduced bispectra 
as a function of scale. In each panel, the solid red lines of increasing thickness denote {Shmmi Bhhmi ^hhh} an d the dashed blue lines 
of increasing thickness denote {Qhmmi Qhhmi Qhhh}- The alphabetical labels {a, b, c. d} correspond to the triangle configurations with 
fci = {0.03, 0.04, 0.05, 0.06} h Mpc -1 , respectively. The vertical black lines denote the effective bias parameter , using the same 

wavemodes that enter into the bispectrum estimates. 



for Bhhm- Hence, we shall think of this as being an almost 
perfect measure independent of discreteness. 

Using these estimators, we compute the ensemble av- 
erage and ensemble-to-ensemble variations of the halo-mass 
cross-bispectra. We do this for the same bispectra configu- 
rations that were considered in £|5 . 1 1 



6.3 Nonlinear bias from cross-bispectra 

We estimate the nonlinear bias parameters and their errors 
from the cross-bispectra using the same method employed 
for the auto-bispectra in H5.2l and H5.3I The results are tab- 
ulated in Tables [5] and [6] respectively. 

Figure \5\ presents the 2-D 95% confidence likelihood 
contours for bi and 62, that are obtained from fitting the 
shot-noise corrected bispectra {-B^hii! Aihmi Bhmm} an d re- 
duced bispectra {Qhhh> Qhhmi QhSm}- The four panels show 
the results obtained from fitting triangle configurations, 
fci <E {0.03, 0.04, 0.05, 0.06} h Mpc -1 , with k 2 /ki = 2, and 
these correspond to the top-left, top-right, bottom left and 
bottom right panels, respectively. For comparative purposes, 
the vertical band in each panel, denotes the lcr constraint 



on b hh ' , obtained from the shot-noise corrected halo and 
nonlinear matter power spectra (c.f. >i4.2|l . 

Consider first the bispectra {-BhSmi -Bff m , BhhLi } 
(solid red lines of increasing thickness), from the figure 
and the tables, we see that all of the results are reason- 
ably consistent with one another over the various scale 
ranges considered. However, when smaller scales are used 
(i.e. fci ^ 0.05 fcMpc -1 ), the consistency weakens and the 
best-fit parameters, obtained from -Bhfii an< A -Bhmmi differ 
by > 2.5cr. 

Evaluating the results for the reduced bispectra 
{QhSmi Qhhmi Qhhh} (dashed blue lines of increasing thick- 
ness), the four panels show a strong evolution of the er- 
ror ellipses with scale. We also note that the level of agree- 
ment between the different estimators also evolves strongly, 
becoming weaker and weaker as smaller scales are consid- 
ered. At the largest scale where k = 0.03 /iMpc -1 , all of the 
2cr likelihood contour regions overlap. However, this consis- 
tency is broken for the next scale range, k — 0.04 h Mpc - , 
where Qhmm and Qhhm are shifted downwards and further 
to the right favoring a more negative &2 and higher 61. The 
trend continues in this same direction heading to smaller 
and smaller scales. 
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Comparing the results from both B and Q together, we 
see that only on the largest scales is there any degree of over- 
all consistency. One way to interpret the results up to now, 
is that, if we believe 61 ~ b^' sc , then the agent driving the 
inconsistency between the parameter estimates, is the break- 
down of the local bias model at tree-level. Furthermore, the 
breakdown of the local tree-level model is more severe for 
the reduced bispectrum than for the bispectrum. 



7 THE NEED FOR BEYOND TREE LEVEL 
BIAS MODELS 

This final set of analysis consists of a simple proof of method 
test, and we determine whether, when the underlying bias 
model is known, the 'true' bias parameters of the model are 
indeed recoverable with our approach. 



7.1 Biasing by hand 

For these tests, and for simplicity, we shall assume that the 
local model of biasing at quadratic order is the correct un- 
derlying bias model. Nonlinear biased density fields of this 
type may be obtained through the following procedure. 

For each of the z — outputs of the 40 simulations, 
we assign the nonlinear density field of matter to a cubical 
Fourier grid using the CIC algorithm. This is then Fourier 
transformed. Each Fourier mode is then smoothed using a 
Gaussian filter of scale R. We then inverse Fourier trans- 
form this field and obtain the smoothed, nonlinear matter 
distribution in real space. Using this we next form the sum, 

St,(x\R) = b\S{x\R) + b b 2 [S(x\R)] 2 /2 , 



(67) 

where b\ and 6 2 are the artificial bias parameters. Finally, 
this is Fourier transformed to give us <5b (k|i?) . Thus, given 
S(k\R) and <5b(k|B), we can now use our standard bispec- 
trum estimators to estimate B b bb, B b b m , B bmm and B mmm . 
We refer to this procedure as the 'biasing- by-hand' test. 

The major benefits of these tests are that we are able to 
better gauge the effects to which nonlinearities beyond tree- 
level order influence the measured bispectra and reduced 
bispectra. We also note that shot-noise plays no role here, 
since the biased field is created from the matter density field 
which is densely sampled. 



where for ease of notation we take b\ = bi and in the above 
we have suppressed the dependence of B, Va,ui, Vs, m and 
Ve,m, on (ki,k2,— ki — k2). We have also introduced the 
auxiliary functions: 



P 4 ,m 



W(k 1 R)W(k 2 R)W(k 3 R)P n , m ■ 



d 3 qi 



IF, 



qi,ki-qi 



(2tt) 3 

xT(qi,ki -qi,k 2 ,k 3 ) + 2cyc ; 

(2tt)3 (2n) 3 cllM ~' 11 ci2M ~ cl2 
xP 5 ,m(qi, ki - qi, q 2 , k 2 - q 2 , k 3 ) 



+2 eye ; 
i 



(72) 



(73) 



(74) 



d 3 qi 



d 3 q 3 



■ W, 



q3> k 3-<13 



(2tt)3 (27r)3" qi ' K1 ~ qi 
xP 6 (qi,ki -qi,q 2 ,k 2 -q 2 ,q 3 ,k 3 - q 3 ) . (75) 

The attractive aspect of this test can now be understood: if 
we move the terms in Eqns (|69|) . (|70|) and <|7ip . which are 
proportional to B mmm from the right to the left-hand-side, 
then we may rewrite this system as the matrix equation: 

6 2 /2 
= | M 2 /3 6|/12 
blb 2 /2 Ml/4 61/8 

where we defined 34mm = Bbmm — 6iB m mm 
tion may be inverted to give, 

6| 



/ 34mm \ 




1 !Vbbm 




V ^bbb ) 







261 







( 3^bmm 


-86i6 2 


12&2 


:)' 


3A>bm 


86? 


-246i 




\ Vhhh 



.(77) 



Hence, if we specify &i, 6 2 and measure the four bispectra 
Bmmm, Bbmmi Bbbm and -Bbbb, then we can determine ex- 
actly 'P4,m, Vb,m and V&,xa- Thus we have complete knowl- 
edge of all components of the nonlinear model at all orders in 
the theory. The lowest order perturbation theory expansions 
of these statistics are (c.f. Appendix |B| : 

6 2 



B 



B 



in) 



&iB 



(0) 



+ 



(0) _ /, 2 r(°) 



-f [v£L(ki)v£L(k 2 ) + 2cyc] ; (78) 
+ ^ [Vn!(ki)V£L(k 2 ) + 2 eye ] ; (79) 
m + bfb 2 [v£L(ki)v£L(k2) + 2 eye] , (80) 
where in the above, we defined Pmm(fc) = W 2 {k\R)P mln (k). 



"bbb ~ °l°n 



7.2 Theoretical interpretation 

In order to interpret the results from such a construction 
we may use the results presented in Appendix [B] with the 
small modification that we do not de-smooth the results. If 
we define the smoothed bispectra as 

S(ki,ka,ks) = W(kiR)W(k 2 R)W(k s R)B(k 1 ,-k 2 ,'k 3 ) , (68) 
then for {Sbmm, Sbbm, Sbbb}, we have: 
6 2 

$b mm — 1 Dm m m H~ 7~^4 . m j 



B h 



Sbbb 



6?B n 



btB 



6162 , 
62. 



6?6 2 



12 

61 &!. 



61. 



m H~ 7~*4,m ~\~ . ^5,m ~t~ Q 7^6, m 5 

Z 4 o 



(69) 
(70) 

(71) 



7.3 Results of the artificial bias test 



Following the algorithm described in t]7.1l for each real- 
ization of our ensemble of simulations, we generate three 
artificially biased density fields smoothed on scales: R = 
{20 , 10 , 6.67} ft _1 Mpc. In all cases we apply the same non- 
linear bias: 61 = 1.63 and 6 2 = —0.53. Whilst these values 
are somewhat arbitrary, they were selected to coincide with 
the best-fit values to the scatter plot of <5h(x|i?) vs. S(x\R), 
smoothed at R ~ 10 /i~ 1 Mpc, that we recorded in i)4.1l 

For each filtering scale, we then measure the four bispec- 
tra -Bmmm, Bbmm, Bbbm and Bbbb for triangle configurations 
with k\ — 0.04/iMpc~ 1 , k 2 /ki — 2, over 20 angular bins. 
From these we use the method described above, to recover 
the higher-order terms: V4,m, "Ps,m and 7\m- 

We now define three modelling cases of interest: 
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R[h- 1 Mpc] 




6i ± 




b 2 ± <Jb 2 


x 2 


on 


-Dbbb 


l.OZ ± 


n n7 
U.U I 


n ap. J— n 1 o 


U.01 


20 


-Bbbm 


1.62 ± 


0.10 


-0.49 ± 0.19 


0.00 


20 


^bmm 


1.63 ± 


0.22 


-0.53 ± 0.51 


0.00 


10 


-Bbbb 


1.62 ± 


0.04 


-0.42 ± 0.04 


0.15 


10 


Bbbm 


1.62 ± 


0.06 


-0.47 ± 0.07 


0.02 


10 


Bbmm 


1.63 ± 


0.13 


-0.53 ± 0.20 


0.00 


6.7 


-Bbbb 


1.59 ± 


0.04 


-0.35 ± 0.02 


0.70 


6.7 


-Bbbm 


1.60 ± 


0.05 


-0.42 ± 0.04 


0.14 


6.7 


^bmm 


1.63 ± 


0.12 


-0.53 ± 0.12 


0.00 



Table 7. Constraints on b\ and b 2 obtained using the Exact 
Trispectrum model described in the text. The actual input 
bias parameters were f>i = 1.63 and b 2 = —0.53. Column 1: the 
smoothing scale of the biased density field; Column 2: measured 
quantity; Column 3 and 4 best-fit values for b\ and b 2 along with 
la errors; Column 5: the median \ 2 ■ 



Rlh^Mpc] 




b 1 ± a bl 


b 2 ± °~b 2 


X 2 


20 


^bbb 


1.63 ± 0.11 


-0.67 ± 0.36 


14.75 


20 


Bbbm 


1.58 ± 0.17 


-0.45 ± 0.65 


15.71 


20 


-Bbmm 


1.37 ± 0.33 


0.56 ± 1.19 


18.43 


10 


Bbbb 


1.49 ± 0.03 


-0.66 ± 0.08 


13.32 


10 


Bbbm 


1.48 ± 0.04 


-0.69 ± 0.13 


13.60 


10 


Bbmm 


1.46 ± 0.09 


-0.69 ± 0.29 


14.13 


6.7 


Bbbb 


1.36 ± 0.02 


-0.74 ± 0.06 


13.21 


6.7 


Bbbm 


1.36 ± 0.03 


-0.82 ± 0.09 


12.87 


6.7 


Bbmm 


1.32 ± 0.07 


-0.86 ± 0.22 


13.13 



Table 8. Same as Table 171 but this time the \ 2 analysis is for 
the Tree-level model described in the text. 



• Case 1: All Order: Eqns (J69J) , (JTQj) and (JTTJ) are used 
to interpret the data. 

• Case 2: Exact Trispectrum: Eqns ([69]), (JTOJ) and (fTI) 

are exact up to 7- , 4, m - All higher-order terms (Ps.m, Ve,m) 
are dropped from the modelling. 

• Case 3: Tree-level: lowest order expansions given by 
Eqns (f78j) . (ff9|) and ((80)) are used to interpret the data. 

For each of the models described above, we then apply the 
same x 2_n tting analysis, as described in H5,2l to determine 
the best-fit bi and 62 parameters. 

We begin by first examining the All Order expansion 
model. We confirm that for this case, the true bias parame- 
ters bi = 1.63 and 62 = —0.53 are recovered exactly, albeit 
with some uncertainty, however, with a \ 2 — 0, and for all 
the smoothing lengths considered. This null test is impor- 
tant, because it gives us confidence that any departures of 
the fits from the true bias values, can be attributed solely 
to a break down of the theoretical modelling. 

Next we focus on the Exact Trispectrum model 
where B mm m and Va, m are measured from the simulations. In 
Table [7J we report the best-fitting bias parameters with the 
la errors expressed in quadrature for the auto- and cross- 
bispectrum and for the four smoothing scales examined. For 
the case Bbmm, a quick inspection of Eq. (|69[) tells us that 
the modelling should be exact, and indeed we see that the 
bias parameters are correctly recovered. However, for the 
cases £>bbm and £>bbb we see that the absence of the higher- 
order terms ("Ps.m, Ve,m), induce biases in the parameters. 



For fei the deviation from the true value is relatively small, 
with the value of the parameter only slightly decreasing in 
size. For 62 the deviations are larger, and this parameter be- 
comes more positive. We also note that the deviations from 
the true values appear to increase as the smoothing scale is 
decreased. 

Finally, we focus on the Tree-level model. Table [8] 
presents the best-fit results for &i and b 2 . We see that in 
nearly all cases, there are systematic biases in the recovery 
of the nonlinear bias parameters for all of the measured bis- 
pectra. In particular, for the case of Bbmm, the results are 
most deviant and poorly constrained. Whereas for Z?bbb, only 
when the data has been smoothed on scales R = 20 /i -1 Mpc 
are the recovered parameters close to the true values. 

The comparison of the results from this analysis leads 
us to conclude that the recovered bias parameters are very 
sensitive to the inclusion of beyond leading order corrections 
in the modelling. Furthermore, accurate nonlinear modelling 
of, at the very least, the matter bispectrum and trispectrum 
will be essential, if we are to safely recover the nonlinear 
bias parameters from this approach. 



8 DISCUSSION 

We have evaluated the local halo bias model at second-order 
using three different probes: smoothed density fields; power 
spectra; and bispectra and reduced bispectra. A summary of 
our results for the best-fitting bias parameters determined 
from shot-noise corrected spectra is shown in Figure [7J 

In the figure we also compare our estimates for bi and 
b 2 with the analytical predictions for the halo bias parame- 
ters obtained from the peak background split (PBS) ansatz 
(|Bardeen et al.lll986l ; IMo fc Whitelll996l ). The average the- 
ory bias parameters are obtained through computing the 
expression: 



* = h 

n 



dMn(M)bi(M) ; n: 



dMn(M) , (81) 



where n(M) is the halo mass function, and Mmm 
is set equal to the value of the minimum halo mass 
identified in the simulations (see $3} . We evaluate the 
above in tegral using three differ e nt fits to N-body sim ula- 
tions bvlSheth fc Tormenl (| 19991 ). IWarren et all (|2006l ) and 
IPillepich et alJ (|2010h . The corresponding expressions for 
th e bias parameters a s a fun ctio n of halo mass are pre sented 
in IScoccimarro et alJ (|200ll ) and iManera et all (|2010l ). 

Considering, the results for bi (bottom panel), we see 
that when the reduced bispectra, Qhhh, Qhhm and Qhmm 
are used, the recovered parameters are poorly constrained 
and appear incompatible with respect to the other estimates 
and are only weakly consistent with one another. On the 
other hand, the estimates from the bispectra Bhhh, Shhm 
and Bhmm are in much better agreement with each other. 
They are also in clo se ag reement with the pre dictions from 
IWarren et all (|2006l ) and IPillepich et al.l |2010l ), which both 
provided an estimate of bi = 1.39. However, they slightly 
undershoot the values from the effective bias estimates, 
fchh' an d &hhi likewise the smoothed density fields, and the 
Sh eth-Tormen prediction. The analytical predictions from 
the lSheth fc Tormenl (|l999l ) mass-function yielded bi = 1.50 
in good agreement with the power spectrum and density field 
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Figure 7. Plot summary of bias measurements on b± and 62 
for the second-order local bias model from different estimators:, 
shot-noise corrected B, Q, P, and finally smoothed, S R , in com- 
parison with analytic al predictions applying th e peak-background 
split ansatz with the ISheth fc Tormer] d 1999T> mass functio n de- 
noted by the dotted line, as well as the lWarren et aT and 
IPillepich et~aH l|20ld l mass-functions, which are both represented 
by the single dashed-line. 



results smoothed on a scale R ~ 50/i -1 Mpc. The recovered 
values of 61 from the effective bias in the power spectra 
and the smoothed density fields collectively are in broad 
agreement, but the latter increase with decreasing smooth- 
ing scale. 

In the top panel of the figure, we see that the con- 
straints on 62 from the different estimators used in the 
simulations are reasonably consistent with one another, al- 
beit with significant error bars. These estimates also agree 
well with the PBS p r edicti on from the I Warren et all (|2006t ) 
and IPillepich et"all (|201dh mass functions, which give an 
average 62 = —0.24. However, the prediction from the 
ISheth fc Tormenl (| 19991 ) mass function gives 62 = —0.06, 
and this appears to be in worse agreement with the data. 

There are a number of possible explanations for the de- 
viations in the recovered bias parameters. Firstly, the rela- 
tion between matter and halo fluctuations may not be local. 
Indeed, we know it is not deterministic owing to the scat- 
ter in relation Sh vs. 5. Perhaps, this is a consequence of 
non-locality. In this case we need a more advanced theoret- 
ical approach to understand the halo clustering. One pos- 
sibility may be that the bias is local in Lagrangian space 
(ICatelan et al. I l2000l ; lMatsubarall201ll ). 

Secondly, our simple biasing-by-hand test has enabled 
us to discern that, for the current set of tests that we have 
performed, the most likely explanation at this point is that 
the tree-level expansions for Bhhh and Qhhh are not suffi- 
ciently accurate enough. Higher-order nonlinear corrections 



in the modelling must be included, and if possible all order 
expansions for -B mm m and P4,m would be invaluable. 

Thirdly, as we have argued, the local bias model only 
makes sense in the context of smoothing. The bias param- 
eters one recovers from fitting, depend sensitively on the 
smoothing scale R. For the biasing-by-hand tests, the exact 
smoothing scale was known beforehand. However, in real 
data we do not know this a priori. In all of the cases, when 
recovering bias parameters from the bispectra, we have as- 
sumed that we are on sufficiently large scales such that 
W(kiR) — » 1. However, in general R should be a free pa- 
ram eter and as such, marg i nalize d over in the analysis. 

anera fc Gaztafiagal (|201ll ) performed a similar study 
of nonlinear halo bias with the three-point correlation func- 
tion in configuration space. In contrast to our analysis, they 
measured the bias parameters for different halo mass bins. 
They found inconsistencies between the predictions of the 
different estimators considered. They evaluated the scatter 
plots of Sh vs. 8 as function of smoothing scale, and found 
that stability in the local bias parameters (61,62) occurred 
for smoothing scales, R > 30—60 h~ 1 Mpc, albeit with larger 
errors. They also found that the bias predictions derived 
from Sh vs. S for R = 60 ft _1 Mpc were in good agreement, 
to within the errors, with the linear bias measured from eval- 
uating the two-point correlation function on large and inter- 
mediate scales. As in the case of our findings, they found the 
linear bias measured from evaluating the three-point corre- 
lation function, expressed in terms of the Q-amplitudes, was 
not consistent with that of the two-point correlation func- 
tion for the lower mass bins: M < 10 13 ^"'Mq. They were 
una ble to formulate sol id conclusions for larger mass ranges. 

iGuo fc Jingl (|2009l ) explored the differences between es- 
timates of bias from Q and P. They also found that 61 based 
on analysis of the mock galaxy catalogues was larger for 
galaxy re duced bispec t ra an d power spectra Q s than for 
P s . While [Guo fc Jin g (2009) also noted that this might be 
due to the failure of SPT at tree-level, they also reported 
that agreement could be found between estimators if Q mmm 
measured directly from the simulations was used in place of 
the tree-level expression. However, when we performed the 
same test with our data we found no dramatic reconcilia- 
tio n of the two b ias es timates. The investigation performed 
by IGuo fc Jingl (|2009h was carried out using only 4 large 
volume and 3 smaller volume runs. As a result of having too 
few realizations, they assumed the Gaussian approximation 
for the covariance matrix in order to perform their study at 
large-scales. 



9 CONCLUSIONS 

In this paper we have used a sample of 40 large volume 
A r -body simulations, with total volume V ~ 135Gpc 3 /i~ 3 , 
to test the local model of halo biasing, and the extent to 
which nonlinearities impact the modelling. We used three 
different methods for exploring the bias: smoothed density 
fields; power spectra; and bispectra and reduced bispectra. 
We focused mainly on the results from the bispectra. All of 
the reported results were scaled to a single realization of our 
simulations, and so are directly relevant for galaxy surveys 
with a total volume V ~ 4 h~ 3 Gpc 3 . 

In ij2]we reviewed the basic results of perturbation the- 
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ory and how they connect to density statistics. We then 
reviewed the local model of halo biasing, drawing special at- 
tention to the role that smoothing plays in the theory. The 
important result being that even at tree-level, the smooth- 
ing explicitly enters the theory. The expressions for Bhhh 
and Qhhh, which are typically used in all past and current 
analysis, make the assumption that smoothing is unimpor- 
tant. In subsequent sections we argued that this assumption 
is not safe. 

In §3]we described our suite of iV-body simulations, and 
the halo catalogues used in this study. 

In ijj] we made measurements of the relation be- 
tween 5h(x|-R) and S(x\R) smoothed on the set of scales 
R = {50, 20, 10}ft _1 Mpc. To this data we fitted the local 
model of halo biasing up to second order, including bo, bi and 
62. We found that the fits were reasonably good, however the 
best-fit parameters showed a running with the filter scale R. 
We then demonstrated, theoretically, why the nonlinear bias 
parameters from this approach, could not be made indepen- 
dent of smoothing scale. We then turned to Fourier space 
statistics, and used the halo auto- and cross-power spec- 
tra to obtain an effective large-scale bias. We found that 
the effective bias estimators b^h and bj^ were reasonably 
scale independent for k < 0.08 ft Mpc -1 . However, on scales 
smaller than this, bj^ decreased with increasing wavenum- 
ber, whereas bh m remained surprisingly flat. 

In iJS] we estimated the matter auto- and halo auto- 
bispectra and reduced bispectra from our simulations. We 
measured these statistics for the triangle configurations 
fci = {0.03, 0.04, 0.05, 0.06}ftMpc -1 and with fc 2 /fci = 2 
and 0i2 £ [0, n] . These triangles all lay in the weakly nonlin- 
ear regime k = 0.03 — 0.18 ft Mpc -1 . We modeled these es- 
timates using tree-level perturbation theory expressions for 
the matter bispectrum and nonlinear bias at second order, 
and assumed smoothing to be unimportant. Our method for 
estimation of the bias parameters followed a standard min- 
imum x 2 approach. We estimated the covariance matrix for 
the full ensemble applying principal component analysis to 
minimize the intrinsic noise. We also performed a jackknifc 
subsampling routine to propagate the error of the estimated 
covariance matrix onto the errors of the bias parameters. 

We tested how well the measurements of the matter 
bispectra B mmm and reduced bispectra Q mmm could be de- 
scribed by such modelling. The results obtained for the bias 
parameters bi and b2 showed that the tree-level expressions 
were a good description of the data for configurations, fci = 
{0.03, 0.04} /iMpc -1 , for which bi = 1 and b 2 = 0. However, 
for smaller scale triangles, ki = {0.05, 0.06} ftMpc" , signif- 
icant deviations were apparent, and these were manifest as 
bi 1 and b 2 ^ at high significance. 

We then applied the \ 2 test to the halo bispectra Bhhh 
and reduced bispectra Qhhh- We found, for the shot-noise 
corrected Bhhh, that the estimated values for bi ~ 1.40 and 
b 2 ~ —0.25 were reasonably consistent with one another. 
However, the fits became progressively poorer as smaller 
scales were added, yet the reduced-x 2 remained < 2 for 
k\ — {0.06} ft Mpc -1 . For the shot-noise corrected Qhhh, 
we found that the values of bi were significantly larger 
bi ~ 1.85, with large errors, and the values evolved with 
triangle configuration scale. However, b2 ~ —0.3, appeared 
to be more stable, although again with large errors. For tri- 
angle configurations with ki ^ 0.04 ft Mpc -1 , the fits from 



Bhhh and Qhhh were inconsistent with each other at the ~ 3<r 
level. For both Bhhh and Qhhh shot-noise corrections signif- 
icantly influenced the recovered bias parameters. 

In ij6]we explored the halo and matter cross-bispectra, 
Bhhm and -Bhmm, and reduced bispectra Qhhm and Qhmm- 
We calculated the tree-level expressions for these quantities 
symmetrized in all of their arguments. We then developed 
estimators for them. We showed that for Bhmm and Qhmm, 
provided the matter distribution was densely sampled, the 
shot-noise corrections were small. 

We applied the \ 2 analysis from jj5]to these statistics 
and recovered the best-fit values for bi and b2- We found 
that for Bhhm the shot-noise corrected data were all rea- 
sonably consistent with one another, giving bi ~ 1.39 and 
b 2 ~ —0.3. For Bhmm we found a similar pattern, except 
that for fci ^ 0.05 ft Mpc -1 where we found b 2 ~ 0.0, but 
with large errors. The results for Qhhm and Qhmm appeared 
to vary significantly. 

Finally in |J7]we explored to what extent the break-down 
could be attributed to the absence of terms that were beyond 
tree-level in the modelling. In order to do this, we developed 
a novel approach, whereby we constructed smoothed biased 
density fields from the smoothed matter density field, using 
the local model at quadratic order. We showed that if we 
set bi and b 2 to some fiducial values, and then measured 
the smoothed matter and halo bispectra and their cross- 
bispectra, then the higher order matter correlators Pi, m , 
P^.m and Bo,m could be recovered exactly. Thus we were 
able to construct three models: an all order model; a model 
that used the exact matter bispectrum and trispectrum; and 
a tree-level model. 

We applied the \ 2 analysis using these three models and 
for bispectra with ki — 0.04 ft Mpc -1 . As expected, the ex- 
act model recovered the correct bias parameters. The model 
with the exact B mmm and Pi, m , was in fact also exact for 
Bbmm- For -Bbbm and -Bbbb the recovered parameters were 
close to the true values, but showed evolution with smooth- 
ing scale. Finally, for the tree-level model we showed that 
there was a significant evolution in the estimated bias pa- 
rameters with smoothing scale and with the type of statistic 
used. 

We conclude that estimates of nonlinear bias from the 
bispectrum that do not attempt to account for higher-order 
corrections, will most likely provide biased estimates for the 
bias parameters bi and b2. Robust modelling of nonlinear 
bias from bispectra, will, at the very least, require almost 
exact models for the matter bispectrum and trispectrum. 

Real space estimates of bias appear to be inconsistent 
with Fourier space based ones. We believe that this owes 
primarily to the mixing of large- and small-scale wavemodes 
in real space. We therefore recommend that perturbative 
methods should strictly be applied in Fourier space. We also 
recommend that measurements focus on the bispectrum and 
the associated cross statistics, rather than the reduced bis- 
pectra, since this appears very sensitive to nonlinearities in 
the modelling and also shot-noise corrections. 

Finally, we emphasize the importance of smoothing in 
the local model. Owing to the fact that the smoothing scale 
associated with the halo/galaxy distribution in question is 
not known a priori, it must be treated as a nuisance param- 
eter and so marginalized over. 

An alternative strategy for recovering information from 
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higher order statistics, which may be of interest for future 
consideration, is the use of 'Gaussianizing transformations ' 
or 'density clipping' (|Nevrinck et al.ll2009l ; TSeo et alj|201ll ; 
ISimpson et aljfeoill ). However, the theoretical connection 
between what is measured and what is interpreted from such 
approaches still remains to be fully calculated. 
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APPENDIX A: BISPECTRUM ESTIMATION: ALGORITHM 



Briefly, the algorithm that we employ is as follows: firstly the dark matter den sity field is computed by ass igning the dark 
matter particles to a cubical grid using the 'Cloud in Cell' (CIC) technique (|Hocknev fc Eastwood! fl988T ). Next the fast 
Fourier transform (FFT) of the gridded density field is computed. Each Fourier mode is then corrected for convolution with 
the Fourier mesh . We do this by dividing out from each mode, the Fourier transform of the window assignment function of 
the CIC scheme (|Hocknev fc Eastwoodlll988l : |jindl2005h : 



sin(7rfci/2fcjv !/ )" 1 2 



7rfci/2fcjv y 



(Al) 



where subscript g denotes gridded quantities, kjss y — nN g /L is the Nyquist frequency of the mesh and iVg is the number of 

Fourier grid cells. 

The estimator for the bispectrum can be written l|Scoccimarro et al]|l998l ): 

where V u is the sample volume (in our case the simulation volume), the normalization factor, Vb, can be written as 
l|Sefusatti et al.ll2006l : Ijoachimi et alj|2009h 

T/ , \ f d^qi d 3 q 2 d 3 q 3 3 D 87r 2 fcifc 2 fc 3 , . , ,3 / aq n 

V B (k 1 ,k a ,ri = J — — — (2.) 5 (q 123 )« (2?r)6 (Afc) , (A3) 

and we write in shorthand J D (qi...n) = <5 D (qi + . . . + q ra ). A practical implementation of the above estimator may be achieved 
through (|Smith et al.ll2008l '): 

B (k 1 ,k 2 ,e 12 ) = ——Jt—— Yl fte[5(k ni )5(k n2 )<5(k_ ni _ n2 )], (A4) 

iVtritfcl, kl, 612) , , 

(ni,n 2 ) 

where superscript d denotes discretized quantities; n; denotes an integer vector from the origin of the fc-space to each mesh 
point; (ni, n 2 ) represents a pair of integer vectors, which lie in thin shells centred on k\ and fc 2 and whose angular separation 
lies in a narrow angular bin centred on #i 2 , and for which k 3 = — ki — k 2 . The upper limit of the sum A^ r j(fci, fc 2 , #i 2 ) 
represents the total number of triangles that have such a configuration. 

The estimator for the bin-averaged reduced bispectrum, Q, is written as: 

- d ^ d ^ dcnom.d 

Q =B IQ , (A5) 

- dcnom.d 

where Q is the estimator for the bin-averaged cyclical terms of the power spectrum generated from first computing the 

bin-averaged power-spectra, Pj . Note that we estimate the power spectra that enter into this product in a slightly different 
way than normal: we use only those modes that go into estimating the particular B triangle configuration to estimate 

_^dciiom,d 

Q (ki, fc 2 , 812). Hence, 

_ d v Nt,i(ki,k 2 ,e 12 ) 

N t vi(ki,k2,0 1 2) , ^ , 

Ol,n 2 ) 

^d 

where i £ {1,2,3} and where P 3 is dependent on the angular bin, since, with |ni| and |n 2 | fixed, we still have 
cos#i 2 = ni • n 2 /]ni||n 2 | and the closure criterion implies n 3 = — ni — n 2 varies as function of #i 2 . Therefore, 

^ dcnom.d ^ d ^ d ^ d ^ d ^ d ^ d 

Q =P 1 P 2 +P 2 P 3 +P 1 P 3 . (A7) 

The estima tes of B d and Q d are then co rrected for discreteness, i.e. shot-noise. For the estimators of interest, the 
corrections are (|Peebledll980l : ISmith et al.ll200Sft : 



shot 



1/n ; (A8) 



B shot = [P 1 +P 2 +P. i ]/n-2/n A ; (A9) 

-'- dcnom d -~- d — d 

Q shot = 2[P 1+ P2+P 3 ]/n-3/n 2 . (A10) 
Shot-noise corrected estimates of the statistics are obtained: 

X = X d - Xshot , (All) 

^- ■ — ^denom ~~ - - dcnom 

where x € jP,B, Q } and where Q = B/Q . Note that the above recipe corrects some typos that are present in 
ISmith et all (|2008r i. 
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APPENDIX B: HALO CROSS-BISPECTRA IN THE LOCAL MODEL 

As was shown in Eq. {TSJ, at quadratic order, the local model of nonlinear biasing can be written as: 
ba(M) f d 3 qi d 3 qi 



8 h (k\R) = b 1 (M)8(k\R) + 



(2tt) 3 (2tt) 3 



<5(q 1 |J?)5(q 2 |7?)(27r)^"(k 1 - qj - q 2 ) 



(Bl) 



where (5(qj|P) = <5(q;)W (q;P) , is the filtered density. Using this model we may now proceed to calculate the halo auto- and 
halo-mass cross-bispectra. 



Bl Halo- mass- mass bispectrum in the local model 

Let us start with the simplest three-point cross-statistic, the halo-matter-matter bispectrum, this can be written: 
(5 h (ki\M : R)5(k 2 \R)5(k 3 \R)) = h(M) <<5(ki|P)£(k 2 |P)<5(k 3 |.R)) 

+ b -^ J §^J^E^ SD ^ " * " <^(qi|-R)5(q2| J R)<5(k 1 | J R)J(k 2 | J R)> . (B2) 

Let us define the smoothed n-point spectra as: 

(5(k 1 \R)...S(k n \R)) = (27r) 3 5(ki + ... + k n )P n (ki,...,k„|i?) , 

= (27r) 3 5(ki + ... + k m )W(fc 1 i?)...W(fc ?l i?)P„(k 1 ,...,k n ) (B3) 

where P 2 = V = W 2 (kR)P, P 3 = & = W(k 1 R)W(k 2 R)W(k a R)B, and where P 4 = f= W(k 1 R)W(k 2 R)W(k a R)W(k 4 R)T. 
We may now integrate over q 2 to obtain 

b 2 (M) r 



(5 h (ki|A/, R)S(k2\R)8(ka\R)) = (27r) d <r (ki + k 2 + k 3 ) 



6i(M)B(ki,k 2 ,k 3 ) + 



J 4^%7"(qi,ki -qi,k 2 ,k 3 , 



(B4) 



2 

On dividing the above expression by W (k l R)W (k 2 R)W (k a R), then we find the halo-mass-mass bispectrum can be written as 



Bhmm(ki, k 2 , k 3 ) = 6i(M)B mmm (ki, k 2 , k 3 ) + b2 ^'^ f ^^■W qi , kl - qi r mmmm (qi, ki - qi,k 2 ,k 3 ) 



2 J (2tt) 3 

We may symmetrize the above result by constructing the sum [Bhmm + Smhm + Pmmh]/3, and this gives us: 



b 2 (M) f d 3 q 



qi, ki - qi, k 2 , k 3 ) + 2cycj 



(B5) 



(B6) 



Bhmm(ki, k 2 , k 3 ) = 6i(M)B mmm (ki,k 2 ,k 3 ) , 

o J (2ti 

On expanding B and T to fourth order in S, the above expression can be approximated as, 

B^ki.ka.ks) « b 1 (M)BSl m (k!,k 2j k 3 ) + MM> [w k2M3 P {0 \k 2 )P^(k 3 ) + 2cyc] . (B7) 

Finally, in the large-scale limit ki — > 0, or for arbitrarily small smoothing scales, kiR <C 1, the above expression becomes, 

b 2 (M) 



B£L(ki,ka,ks) « 6 1 (M)B^ m (k 1 ,k 2) k 3 ) + [p^(fca)P^(*s) + 2cyc 

The reduced bispectrum Qhmm is given by 
Bhmm(kl, k 2 , k 3 ) 



PPhn 



<5hmm(ki,k 2 ,k 3 ) = 

where 

2 1 

PPhmm = - [Phm(fcl)Pmm(fc 2 ) + 2 CyC ] + - [P hm (fcl )P hm (fe 2 ) + 2 CyC ] 



(B8) 



(B9) 



(BIO) 



In order to calculate the reduced halo-mass cross-bispectrum, then we need to evaluate the halo-matter power spectrum. 
In the local model and up to quadratic order in the bias we have, 

62 (M) f d 3 qi 



PhmO) = 6i(M)P mm (ki) + 



2 J (2tt) 3 
Using the above expression in Eq. (|B 10|) we find 

2 " 



- 3 -W / qi,ki-qiPmmm(qi, ki — qi, — ki) 



(Bll) 



PPhmm(kl,k 2 ,k 3 ) = 



6i(M)P mm (ki) + f 4%W qi , kl - qi B mmm (qi,ki - qi, -krj 



1 

+ 3 



6i(M)P mm (ki) + 



2 J (2tt) 3 
b 2 (M) /• d 3 q 



Pmm(fo) + 2 CyC 



2 y (2tt) 

|3 



3 VKqj^.k! - qi Pmmin(qi, kl — qi, — kl 

+ 2 eye 



^2 /■ d qi 

6i(M)P mm (k 2 ) + — J ■^p-W / q2 ,k 2 -q 2 Bmmm(q2,k 2 - q 2 , -k : 

If we expand PPhmm to fourth order in the density then the above expression simplifies to: 
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^hnL(ki,k 2 ,k 3 ) « [2 + 61 (M)] [pi^CfcOpW (fe) + 2cyc] 



Hence we have, 



<L(ki,k 2 ,k 3 ) 



2 + &i(M) 



Qmmm(kl,k 2 ,k 3 



6 2 (M) 



2&i(M) + &?(M) 



W^^PfitblJffilb) + 2 eye 



10) 



P4 t ii(fci)Pi t i(fc 2 ) + 2cyc 



Finally, in the limit that kiR <C 1, W —¥ 1 and the above result can be approximated by 



Q£L(ki,k a ,k3 



2 + 6i(M) 



Qmmm (kl , k 2 , k 3 



62 (M) 



26i(M) +6?(M) 



(B13) 



(B14) 



(B15) 



B2 Halo-halo-mass bispectrum in the local model 

Again, using Eq. (|B1[) . the halo-halo-mass bispectrum, symmetrized in the fcj arguments, can be written: 



-Bhhmfkl, k 2 , k 3 



b 1 (M)B ki,k 2 ,k 3 + 



&i(M)6 2 (M) /• d 3 q 



[w qi ,k 1 - qi T , mmmm(qi, ki - qi,k 2 ,k 3 ) + 2 eye 



+ - 



&!(M) /■ d 3 qi d 3 q 2 



3 7 (2tt) 

[Wqi,k 1 -q I W q2 ,k 2 -q 2 -P5,m(qi,ki - qi,q 2 ,k 2 - q 2 ,k 3 ) + 2 eye 

th 



(B16) 



12 J (2tt) 3 (2tt) 3 

If we use perturbation theory to expand P, B, T and P5, and only keep terms that are 4 tn order in the density field, then the 
above expression can be approximated by: 

fl^ki.ka.ks) « 6?(M) J B^ m (k lj k 2j k 3 ) + i&i(M)ba(Af) [w kl , k2 Pil(fci)Pi°i(fc 2 ) +2cyc] . (B17) 
In the large-scale limit kiR — > 0, the above expression can be approximated as: 



P^ m (ki,k 2 ,k 3 ) » 6f(M)S^ m (k 1 ,k 2 ,k 3 ) + -6r(M)6 2 (M) |PiZ(M^(fc 2 ) + 2cyc 

The reduced bispectrum is given by, 
Phhm(ki, k 2 , k 3 ) 



PPhh 



Qhhm(ki,k 2 ,k 3 ) = 
where 

PPhhm = | [Phh(fci)Phm(fc 2 ) + 2cyc] + i [P hm (fci)P hm (fc 2 ) + 2cyc] . 

The halo-mass power spectrum is given by Eq. (|B11[) and the halo auto-power spectrum is given by: 



Phh(fc) = 6?(M)P mm (ki) + 6i(M)6 2 (M) 



d 3 qi 
(2ny 



H / q 1 ,k 1 -q 1 Pmmm(qi, kl - qi, -kl) 



6|(M) f d 3 qi d 3 q 2 ~ ~ , . , 

H — J (27r) 3 (27r) 3 qi- k i~qi ^q2,-ki-q 2 -'mmmm(qi, ki - qi , q 3 , -ki - q 2 J 

Expanding the above expression to lowest order in perturbation theory gives, 



Z - 2J ^ 1 [26i(M) + 1] [p£l{k,)P£l{k,) + 2cye 



PP, 

Using the above expression, we find that the tree-level expression for the reduced bispectrum can be written: 



QSL(ki,k 2 ,k 3 



2&i(M) + 1 



Qmmm(kl,k2,k 3 ) + 



2b 2 (M) 



26? (M) + bi(M) 



W kl , k2 P£l(ki)P£Uk2) + 2 eye 



(0) . 



P4° I l(fci)Pi u I i(fc 2 ) + 2cyc 



;(<» 



In the large-scale limit kiR — > 0, we again have W — > 1 and 



QSLCki.ka.ks 



26i (M) +1 



Qmmm(kl,k2,k 3 ) 



26 2 (M) 



26f(M) + 6i(M) 



(B18) 
(B19) 
(B20) 



(B21) 



(B22) 



(B23) 



(B24) 



B3 Halo- halo- halo bispectrum in the local model 

Again using Eq. (|B1|I . the halo-halo-halo bispectrum, symmetrized in the ki arguments, can be written: 



Phhh(kl, k2, k 3 



hi 



+-blb 

4 



(M)B mmm (ki,k 2 ,k 3 ) + h,\{M)b 2 {M) j ^-51 [w qi , kl _ qi T mmmm (q 1 ,k 1 - qi, k 2 , k 3 ) + 2 eye j 
W qi , kl _ qi Wq 2 , k2 _ q2 P5,m(qi, ki - qi , q 2 , k 2 - q 2 , k 3 ) + 2cycj 
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, 2 f d 3 qi d 3 q 2 
2 ' (2tt) 3 (2tt) 3 
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+ 1T 



6| /" d 3 qi d 3 q 2 d 3 q 3 



(2tt) 3 (2tt) 3 (2tt) 3 



W / qi,k 1 -qi^q2,k2-q2W ; q3,k3-q3^ , 6,m(qi,ki - qi,q 2 ,k 2 - q 2 ,q3,k 3 - qs) 



(B25) 



If we use perturbation theory to expand P, B, T, P5, and Pe, and keep only terms that are fourth order in the density field, 
then the above expression can be approximated by: 



B^(ki,ka,ks) ~ &5(M)B<2L(ki,ka 1 ka) + bl(M)b 2 (M) [w* lM P<£l(ki)P£l(k2) + 2cyc] 

In the large-scale limit kiR — > 0, the above expression can be approximated as: 

S^(k 1 ,k 2 ,k 3 )«6?(M)S^(ki,k 2 ,k 3 ) + 6f(M)6 2 (M) [i^ (Ai)P^(*a) + 2cyc] . 

The reduced halo-halo-halo bispectrum is given by, 
Bhhh(ki,k 2 ,k 3 ) 



PP h 



Qhhh(ki,k 2 ,k 3 
where 

PPhhh = [Phh(ki)P hh (k 2 ) + 2 eye 



(B26) 
(B27) 

(B28) 
(B29) 



where the halo auto-power spectrum is given by Eq. (|B21[) . On expanding the above expression to fourth order in the density 
we find, 



PP 



(0) 



bt(M) [pW (fci)pW (fa) + 2 eye] 



Using the above expression, we find that the tree-level expression for the reduced bispectrum can be written: 



Q<Z(ki,k 2 ,k 3 



-L- o<°> rki k 2 ko + 



Wx 1>k2 P£L(ki)P£l(k2) + 2 eye 



,(0) 



(fci)p£i(fc 2 ) + 2cyc 



In the large-scale limit kiR — > 0, we again have W — ► 1 and 



1 



61 (M) 



Qmmm(kl,k 2 ,k 3 ) + 



b2(M) 



(B30) 



(B31) 



(B32) 
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